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Fully Plastic Rotating Disc with Large Strains’ 


MELVIN ZAID?{ 
Sperry Gyroscope Companyt 


ABSTRACT 


The thin rotating disc of arbitrary profile and boundary con- 
ditions is considered for a material with an arbitrary stress- 
strain curve. In order to investigate stresses, strains, and speeds 
leading to failure, a finite strain theory is used. 

The concept of a strain plane is introduced, and with it, general 
solutions are obtained and particularized to many specific pro- 
files and materials. Residual stresses in one particular case are 
evaluated. With the aid of these solutions, it is shown that the 
bursting speed is dependent largely on profile shape, and examples 
are presented to show this relationship. 

The small strain theory indicates that complete yielding of a 
perfectly plastic disc takes place for a speed of about 17 per cent 
above that for a disc at the elastic limit. The large strain theory, 
however, indicates that a speed of approximately 100 per cent 
above the elastic limit can be achieved before the disc will burst. 
The latter is a more realistic theory for evaluation of allowable 
safety factors. 

A theoretical solution for bursting is checked against previously 
performed experimental work and agrees to within 1 per cent. 
Strain distributions check quite well. 


NOMENCLATURE 


instantaneous radius to edge of disc 
original radius to edge of disc 
“DJ — (dr — * | 


Di 1 — o + Cc — — 


undetermined radius 
D(2 — ¢,|** 
constant 
c/b 
variable coefficient relating stress and strain ** 
modulus of elasticity 
es (In ty)** 
2 D dx» 
= deformed length 
= undeformed length 


| ee 
Tx, + 


= to at Xp 


Received January 9, 1953. 

* Based on Doctoral Thesis, Massachusetts Institute of Tech- 
nology, 1951. The author is deeply indebted to Dr. C. W 
MacGregor for his invaluable aid while at M.I.T. 

t Senior Engineer. 

t Division of The Sperry Corporation. 

** A further subscript 0, i, or 1 refers to the outside or inside or 


Position only slightly removed from 7. 


radial distance to an element in the undeformed 
disc 
r/a** 
rja’?? 
r/ b** 
effective stress = 
(V2/2)V(s1 — 52) + (se — 53)® + (53 — 51)? 
effective stress at elastic limit 
principal stresses 
radial, tangential, and axial stresses, respec- 
tively, corresponding to the position (r + u) 
in the deformed disc 
deformed thickness 
undeformed thicknesstt 
w*b*y/So 
w*a?y/So 
radial displacement of an element at initial ra 
dius r during plastic yielding of the disc 
r + u** 
x/a** 
x/b** 
density (mass per unit*) 
fO,n #k 
lLa=z=k 
effective strain = se 
(V2/3)V (4 — e)? + (e — 6)? +(e — a)? 
radial, tangential and axial strain, respectively, 
corresponding to an element at distance (7 + 
u) in the deformed disc** 
angle 
normalizing factor 
s/So 
Sr/So, S¢/So, Sa/So 
oc,’ Cresidual 
op = dOplastic 
Ge = Gelastic 
Trz, Trt) Tee = Shear stresses 
or 2e, + «:** 
di 2e: + " 
w angular speed (rad. per unit time) 


wo = bursting speed of flat disc 


tf A second subscript ¢ refers to the thickness at the edge. 
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(I) INTRODUCTION 


HE THIN ROTATING DISC is often met in practical 
"hl ace problems. If the speed is sufficiently low, 
the disc will behave elastically, in which case the solu- 
tion is known.' For intermediate speeds that lead to 
consideration of the plastic-elastic disc, the solution is 
also known.?* However, except for experimental 
work, little has been done on a large strain theory for 
the thin rotating disc. Such a theory is useful in pre- 
dicting bursting speeds and thus in establishing safety 
factors. 


Since this paper represents an early step toward the 
development of such a theory, it is felt that certain 
For instance, a 
From 


simplying assumptions are justified. 
deformation rather than a flow theory is used. 
a purely mathematical point of view, it is sufficient to 
state this assumption and then solve the resulting equa- 
tions. However, such a procedure can lead to unrealis- 
tic solutions. It is therefore worth while to state in 
advance some of the results obtained for this problem. 
(1) The only case for which a detailed study is made 
shows that the stress ratios remain nearly constant as 
the loading is increased. (2) For another case, in 
which the results had been previously established ex- 
perimentally, the bursting speed is predicted to within 
1 per cent, and the strain distributions agree quite well 
even for strains beyond the maximum load point of the 
tensile specimen. 

Solutions are obtained for discs with flat and variable 
profiles and for a strain hardening or nonstrain-harden- 


ing material. 


(II) THEORY 


It is assumed that: 

(1) The material is homogeneous and isotropic. 

(2) The sum of the plastic strains equals zero. 

(3) The ratio of the principal true strain differences 
to the corresponding stress differences are equal. 

(4) There is a unique relationship between effective 
stress and effective strain. 

(5) The variation of stress across the thickness is 
neglected, as well as the effect of axial components of 
the radial stress for large deformations. 

(6) Heating due to air friction is negligible. 


The effective stress at a point is defined as 
S = (VW2/2)V (5, — 52)? + (52 — 53)? + (53 — 1)? (1) 
and the effective strain, as 
e= (v/2/3)V (a — a)? + (2 — a)? + (a — 4)? (2) 


Assumption (4) can be derived from a tension test in 
the form 


s = f(e) (3) 
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For the rotating disc 
5; = S; Gg = & 
Sg = Sy €2 = €& 
S3 = 5. = 0 63 = € 
and we take 
Tre = Te Tz = 0 


Assumptions (2) and (3) lead to 


Ee, = D*[s, — (1/2)s;] (4 
Ee, = D*[s, — (1/2)s,] (5 
Ee, = —(D*/2)(s, + 5,) (6 


The true strain is defined, after Ludwik,’ as 


de = Gi/i 
or 
e = In (L/h) 
This leads to 
I(r u) 1u 
a =in(1 +“) (7 
dr dr 
(r + u) dé u 
e = in = ai In (1 + ) (8 
r dé r 
& = In (t/to) (9) 
The equation of equilibrium taken in the radial direc- 
tion is 
ds, St — S, d 
; = — S (In t) — w?y(r + u) 
d(r + u) r+u d(r + nu) 
(10) 


Combining Eqs. (7) and (8), one obtains the compatibil- 
ity equation, 


1 — ef ~* = (r+ u)[de,/d(r + u)] (11) 


It is now convenient to introduce the following di- 
mensionless variables: 


D = D* se, 


c= 


xX = (r + u)/b 
Se) Si, o, = S,/So 
T = w*b*y/s0 

Where + is an undetermined radius and so an unde- 
termined stress, usually taken as the elastic limit. 
The equations then become 


= = = z ee O; ™ (In t) — Tx, (12 
Xp(de,/d%) = 1—e*° * (13 

o = g(e) (14) 

207 = (6, — o,)*?* + 6? + o? (15 
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FULLY PLASTIC ROTATING 
e. = Dio; — (e,/2)] (18) 
e, = —(D/2)(o, + o;) (19) 
e, = In (t/to) (20) 
Let us define two new symbols 
go, = 26, + & (21) 
od; = 2e, + « (22) 
where ¢, — ¢, = €: — &. Then Eqs. (17) and (18) 


can be combined to yield 


a, = (2/3)(¢,/D) (23) 
o, = (2/3)(¢,/D) (24) 

Using these in Eqs. (15) and (16), there results 
e = Do (25) 


Rewriting Eq. (12) with these new quantities results in 


¢,D' + D k at op eee 0’ | — 
Xp 
on 
< D*Tx = 0 (26) 
where 
(d/dx,)( )=()/) 


Eqs. (13) to (20) inclusive and Eq. (26) supply nine 
independent equations with eleven unknowns, but, by 
assuning fo, these are reduced to a system that can be 
solved in terms of one variable, x,. It should be noted 
that Eq. (25) is derived from the other nine equations 
and therefore is not independent. 


(III) SoLuTION 


The method of solution will be to write this group as a 
simple pair of linear equations that can be solved nu- 
merically. 

Introducing Eqs. (13) and (20) into Eq. (26) there 
follows: 


de, 


gt or 


3 D?Tx,? — D(¢d, — ¢,) + 


l—e 


1 
Do,X» (In | = ¢,dD — D dd, — Dd, de, (27) 
aX», 


where we used 


d 


- (In fy) 
dx, . 


de, 
(in?) = 4 
dx» 


dX» 
Using assumption (2), there results 


A de, + Bde, 


(28) 


é,adD = 


where 
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Fic. 1. Boundary conditions on strain plane 
Xx DJ ia (op = ?,) 
A = Df -#+- gong 
l - e G 

B = D[2 — ¢@,] (29) 

S o; 
J=-Im+ (In to) 

oC  -_- 


To get the second equation, we write Eq. (16) in its 


differential form, 


(3/2)ede = 6, de, + o, de, (30) 
Eqs. (28) and (30) are now written as 
(3/2) Bede — ,*dD 
de, = 
Bo, — Ag, (31) 
de, = (¢,dD — A de,)/B 


Eq. (31) establishes the strain 


(Bd, — Ad, ~ 0). 
path on a strain plane (Fig. 1) prescribed by Eq. (16). 

For the special case of Bd, — Ag, = 0, it can be 
shown that the strain path is tangent to the lines of e = 
constant or 


—¢,/o; = —B/A (32) 


de.Jde = 


A second special case occurs when the strain path 


crosses the line e, = ¢,, and x, ~ 0. It can be shown 


in this case that 








372 JOURNAL OF 
(3/2)Be de — o,2 dD 
‘em — x,D*J@, (e, — €,) (33) 
or de, > 0, (€, > €,). Using Eq. (33) in Eq. (13) 
(3/2)Be de — $,? dD : 
dX» —> —" = - ¢,D2J (34) 
de, > (¢,/B) dD (35) 


The procedure will then be to use Eqs. (13), (25), 
and (31), or any of the equations for special cases, to- 
gether with the stress-strain diagram, to numerically 
evaluate the strain and instantaneous position (x). 
The original position of the point in question is then 
found from Eq. (8), 

f= HSC (36) 
It will thus only be necessary to establish the boundary 
conditions in strain components. Then, once a set of 
starting values is determined in the strain plane, a path 
is obtained having its terminal point at an intersection 
with the boundary line. The sign of + |de| is selected 
so as to keep dx, single valued, a physical necessity. 


(IV) BoUNDARY AND STARTING CONDITIONS 


Because the numerical work is made much easier, all 
integration will be started at the center of the disc and 
will progress outward. 

(1) Solid Center 


From the equilibrium Eq. (12) we see that 0, — o, = 


0 in order that stresses remain finite as x, ~ 0. This 
means in Eqs. (23) and (24) that ¢, = ¢,. 
Using this in Eq. (16) gives 
e = 2e, = 2e, (37) 


at x, > 0. 


To determine the starting increment of strain on the 
strain plane, we return to Eq. (13). Evaluating the 
indeterminate at x, — 0, making use of the fact that 
€- = €,, leads to 


Ide; = de; (38) 

and using this in Eq. (30) finally results in 
de, = (1/8)de, de, = (2/3)de (39) 
at x» ~ 0. The evaluation of dx, is obtained by dif- 
ferencing from a point x», toward x, = 0. In this 
backward difference we use the values of da, = —de, 
de, = —de,, dey = —de,, and dD, = —dD, the values 


of dD, de being determined from the o — ¢ diagram and 
de,, de, from Eq. (39). Using Eq. (28) leads to 
3€1(dD,/den) = A, + 2B, 


Substituting the values for A;, B, from Eq. (29), con- 


THE AERONAUTICAL SCIENCES 


JUNE, 1953 


sidering ¢; ~ @n, and applying the compatibility 
equation to evaluate x»1/(¢@n — $n), we have 
XD); = —3B, de; + or dD (40 
Egs. (37) to (40) furnish all the necessary starting coy 
ditions for the solid center. 
(2) Free Edge 
Here we have the boundary condition of ¢, = 0 or 
¢, = 0 (41 
x» is chosen in an arbitrary fashion at the inner edge 
By use of a normalizing procedure to be developed, w 
will see that this choice results in no loss in generality 
(3) Rigid Edge 


This condition implies that the boundary suffers n 
change of radius or 


e, = 0 (42 
As above, an arbitrary choice of x, is made. 


(4) General Edge Condition 


Any condition of stress or strain may be imposed at 
the boundaries and, with the aid of Eqs. (23) and (24 
may be expressed in terms of strain alone. 


(V) STRAIN PLANE 


Fig. | illustrates the strain plane, with the lines oi 
constant the abscissa 
€,/€:i, and the ordinate being e,/ €;;. 
the coordinates of the order of 1, the divisor e€,; was 
used, this value being the tangential strain at the center 
of the disc. In the case of rigid center we use e,; instead 
of «,, which is zero. 


e/e,, being elliptical, being 


In order to kee; 


(VI) NORMALIZATION OF T AND x» 


x» is defined as (r + u)/b, where b is undetermined 
Let us further define 


(r + u)/a, Xa = 1 (43 


Xq = 


where a is merely the instantaneous outer radius. We 
then have 


Xa = Xo/Xo0 (44 


as the normalized instantaneous radius. To get the 





original radius, we use Eq. (36), which is written in 
terms of the normalized instantaneous radius as 


Ya = Xa e (45 


where 7, = r/a. 
It will be seen from Eqs. (12), (29), and (31) that 

this normalization requires a change of T to 
T; = 


wa*y/so 


or 
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factor 1 
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PLASTIC ROTATING 


FULLY 


7; = Xno7l (46) 


where we call 7, the characteristic speed of the nor- 


malized disc. 


(VII) GENERALIZED PROFILE 


Several methods* for analysis of the arbitrary profile 


have been discussed. The objective here is to evaluate 


(d/dx») (In ty) 


Using Eq. (7), we have 
dx,/dr, = e 


Therefore, 


d a 
(In ty) (47) 


(In tf) = e 
dry 


AX»y 

For any special profile in which there appears fy = 
fir), Eq. (47) can be used directly. 

In this respect a profile can be constructed, in a piece- 

wise manner, either graphically or analytically. 


A useful profile is 


ik = const., 0<% <n (48) 
i = e's, 11 <% <2 (49) 
fo = const., %2< 1% < fo (50) 


then, 


0 < % <K; 


(in fg) = O, 
UXy 


€ 


d 
(In fo) = ne “, 
Xp 


Xx, Me Ke Ve 


ad 
(In fo) = 0, 72 
aX» 


< % < Te 


which is a decreasing or increasing profile, depending 
on, with a flat at the center and edge. 
In this paper a second profile was studied. Let us 


take the most general case 


to = pg{1 — 2 A,6"] 
= 
where 
B = (ry — Tp1')/(Cyo — Toi’ 
with c, = c/b > 1, 7,’ a constant. 
A useful simplification is to take A, = 6,”, so that 


lo ie p(l — B”") (51) 


Here, g is the thickness of the disc at x, = x, and pisa 
factor required to normalize f)/g with reference to this 


position. From Eq. (51), 


, 
Cit, ~— Tei J 


(> = my 
= / 
i —- 


i 


(In t)) = 
dr, 
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ae 
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Introducing this into Eq. (29) and using Eq. (18), we 


have 


(In f) = - (52) 
AX»p Xvi 
oii evi’ 
le / 
Vp Xpi 
e evi 


From Eq. (52) we see that the factors pg do not appear 
in the solution, thus this solution is general with regard 
to thickness. Hence, it is sufficient to find the solution 
for a normalized thickness and multiply it by g for the 
actual disc. 

To find the normalizing factor p, we need only set 


to/q = lat x, = Xi, OF 
Y> Xi 
et etti’ 
ee — (53) 
Xvi 
~~ ett’ 


e“ is taken arbitrarily, and in this paper 
< ry; with the 


Now fa:’ = Xox’ 
the only case investigated is 0 < 7)’ 
special restriction of rp;’ < ri, for n < 1, in order to 
avoid infinite derivatives. 
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to/q = 0, so that c is taken to be the measurement from 
the center of the disc to the vertex of the extended pro- 
file. This is illustrated in Figs. 2 and 3. 

Now, for ¢ = ©, we have the flat profile. 
1, a linear variation of thickness with radius results. 
Then, since 8 < 1, we have 1 — 6” S 1 — 6,n S 1. 
With respect to the linear profile (7 = 1), we see that 
n < 1 leads to a concave profile and 1 > 1 to a convex 
illustrates these profiles. 


For 2 = 


Fig. 2 


profile. 


curve is shown in reference 8. For purposes of evalua- 
tion it is well to calculate D = f(e) by use of Eq. (29), 
numerically smoothing this data to avoid erratic be- 
havior in the difference equations. 

Values of €; and 7 are then assumed, and the solution 
is obtained numerically with the aid of Eqs. (13) and 
(31). T = 10 is a convenient assumption. The us 
of a flat profile leads to d/dx, (In to) = 0 and the result- 


ant simplification of J in Eq. (29). The starting cond- 
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tions are given by Eqs. (39) and (40), and the final 
condition of ¢, = 0 is shown in Fig. 1. 

Fig. 4 illustrates a typical strain path, while Figs. 
5, 6, and 7 show strain, stress, and half-profile distor 
tion the instantaneous radius. 
More results are shown in the original thesis.‘ 

During the course of the investigation it was found 


versus normalized 


that, if a small error were made at one point in the 
numerical evaluation the next line would compensate 
almost completely so as to wipe out this error. This 
happened so consistently that it is felt the equations 
must lead to a stable solution, which is important in 
numerical work. To test this idea, a large error was 
introduced early in the calculations, and the work was 
then correctly completed. Although the strain error 
did not completely vanish, the relative error in x, did 
sensibly vanish. Since the most important quantity 
is 1}, which is only dependent on x,, we can conclude 
that even a reasonably large error in calculations should 
not affect the final result adversely. 

Fig. Sis a plot of 7; (normalized speed) versus the 
strain at the center of the disc (¢;). Unrestrained plas 
tic flow is taken to be the maximum value of 7; or to 
occur when ¢€; increases without an increase in speed. 
In the case considered, 7, ~ 9. 

Fig. 9 shows a plot of o,/c, vs. ¢; for elements at 


Various initial radii. To derive this, it is first necessary 
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10. w/w = 0.1916, e, = 0.5. 


FIG 

5 pice ‘ ; 
to plot rq: vs. Xq, Where a’ is the initial outside radius 
of the disc. Then for various values of r,, = constant, 
X, is obtained and, hence, the stress ratio at the par- 


ticular e;. To evaluate r,’, write 7ge(a/a’) = 1. Thatis 
at the rim 77 = a’, where 7,9 = 7%o/a. Then, 

a/a’ = (1/xa0)e = e%*, 
since Xz = 1. But 

r,(a/a’) = fq = X%q(e*/e*) (54) 


Taking various values of 7,, = constant for different 
e, (different loadings), we see in Fig. 9 that o,/o, ~ 
constant and oa,/0o, 0. Thus we expect the flow and 
deformation theories to yield nearly the same results. 


(ii) Variable profile: In the variable profile and cen- 
tral hole solutions, it was decided to minimize the work 
by solving only for one ¢;. This allows for a more 
concise investigation of various profiles and inside to 
outside diameter ratios. In value of 
sufficiently large to ensure a bursting speed was taken. 
In an actual case this would finally be checked by tak- 
ing a slightly larger and smaller ¢; to confirm this fact. 
Reference to Fig. 8 will show that this means e; must 
constant, independent 


all cases a €; 


be large enough to make 7) ~ 
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of e;. With this in mind, e; = 0.5 was selected for all 
solutions. 


Fig. 10 indicates the results as w/wp vs. foe/q, where 
w is the ‘“‘bursting’’ speed of the variable profile, w, is 
that of the flat disc, f, is the rim thickness, and g is 
the center thickness. The curve applies to the linear 
case with two points for the concave (m = 1/2) and 
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convex (7 = 2) profiles. It is seen that the removal oj 
metal closer to the rim increases the bursting speed. 

(B) Steel._-Published in a paper by MacGregor and 
Tierney" are results of experimental work with high. 
speed rotating discs. Only flat discs were investigated, 
and one of them, identified as disc No. 11, had essen. 
tially the boundary conditions of solid center and stress. 
free outer edge. The material was steel, of which the 
average true stress-strain curve of many specimens js 
shown in Fig. 35 of reference 8. The disc diameter was 
26 in., and the bursting speed was 11,800 r.p.m. Strain 
patterns of that paper are reproduced here for the pur- 
pose of comparison. 

Fig. 11 shows a curve of w X 107° (r.p.m.) vs. €; with 
the speed for unrestrained flow being 11,700 r.p.m. 
This compares within 1 per cent of the experimental 
figure obtained. Fig. 12 compares the experimental 
and theoretical strain distribution. Other comparisons 
are given in reference 8. 


(2) Center and Edge Stress Free 


(A) Brass.—(1) Flat disc and varied profile: Fig, 
13 illustrates the graph of 7) vs. rim to hole ratio. 
The x;¢/x; is the instantaneous value, whereas 7/7; is 
in terms of the original radii. This graph includes the 
flat and variable profile. As with the discs lacking a 
central hole, the removal of metal increases the burst- 
ing speed. Also, there seems to be no stress concen- 
trating effect, but instead the flat disc with central 
hole shows a definite tendency toward the solution for 
the disc without central hole as the hole diameter de- 
creases. 

Fig. 14 shows typical strain paths. 

(11) Thin ring: From Fig. 14, it is seen that, as the 
ring thickness decreases, the strain path decreases and 
becomes tangent to ¢, = 0. This is also apparent from 
the fact that for the thin ring ¢, = constant throughout. 
Taken with the boundary conditions, ¢, = 0 every- 
where in the ring. Further, as in previous work, the 
integration is started at ¢, = 0 and stopped as soon as 
it again reaches ¢, = 0. Thus the strain path lies on 
the line ¢, = O and is specified by two neighboring 
points that approach each other in the limit. 


From the above and Eq. (28), 
de, = (—A/Bhde, 
2e, = —e, (everywhere) 


Taking the differential of the second equation, we 
find 


A/B = 1/2 
Introducing this into Eq. (29), there finally results 
Tx,? = ¢/D = ¢6/D=c¢ 
But, from Eq. (46) and the fact that x,, = x = *w 


T;, = w*a*y/s5 = o (90 
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or 


wa*y = S$ (56) 


where a is the instantaneous radius of the ring. For 
small strains, a > r and Eq. (56) reduces to the known 
equation for a thin ring.® 

In order to find the original radius r, it is only neces- 
sary to use Eq. (36), which can be written as 


€ 


a= re (57) 
while the change in thickness becomes 
= » - 
t/tt =e *° (58) 


irom Eq. (20) and assumption (2). 

Eqs. (55) to (58) inclusive completely specify the 
free thin ring. 
t/tt = 0.779 and 7; = o. 


As an example for € = 0.5, x/r = 1.649, 


(tii) Ideally plastic material: For this « = Day, 
Thus, de = dD, which 
Examples have 


where o9 = constant = 1. 
slightly simplifies the calculations. 
been worked for this case and are shown in detail in 
relerence 8. Comparison with the small strain solu- 
tion of Nadai and Donnell? shows similar stress and 


strain diagrams. 


(3) Residual Stresses 


The problem of finding residual stresses is usually 
complicated; however, if unloading and reloading is 
assumed to be perfectly elastic without a hysteresis 
loop and if the Bauschinger effect is ignored, these 
stresses can be evaluated. In this case, the residual 
stresses are equal to the plastic stress less the stress 
obtained by assuming elastic behavior of the disc at 
the same speed. 


Cy! = Gy — Ce (59) 


Only one such case was investigated, and this was 
taken for a disc that becomes flat at the prescribed 
speed. The elastic solution for this case is then greatly 
simplified and well known.' 

For the plastic solution, ¢ = constant will lead to 
further simplifications in the equilibrium equation and, 
hence, in the basic Eq. (31). These are basically in the 
coefficient of A and B. Such an analysis is shown in 
Fig. 15 illustrates the result of a 
central 


the original thesis.* 


calculation for these stresses for brass with 


Strain e; = 0.22. 
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An Investigation of the Flow About Cones 
and Wedges at and Beyond the 
Critical Angle 


G. W. JOHNSTON * 
Institute of Aerophysics, Unwersity of Toronto 


SUMMARY 


An experimental investigation has been carried out in the 
Institute of Aerophysics’ 16- by 16-in. supersonic wind tunnel 
to determine the extent and general nature of the subsonic field 
produced ahead of wedge and cone models during bow-wave 
detachment and beyond, where the free-stream Mach Number 
lies in the supersonic régime. All tests were carried out at a 
nominal Mach Number of 2.5 lift 
The extent and growth, both lateral and forward, of the im 


under conditions of zero 
bedded subsonic field has been measured for both models at the 
critical angle and beyond from a series of schiieren and shadow 
flow photographs. The general nature of the associated sub 
sonic field was determined over the same range from a series 
of pressure distributions taken at the model surface and from 
the measured The 
ments have also been used to determine the forebody pressure 


shock-wave geometry. pressure measure 
drag of cones and wedges at the detachment angle and through 
out the critical position of the detached region. A method has 
been developed which allows a complete theoretical solution for 
the detached subsonic flow about a wedge model This method 
rely on a preliminary knowledge of the bow-wave 
geometry but is carried out ii ihe hodograph plane. 


The salient characteristics of the detached flow about 


does not 


both 
the cone and wedge models are shown in Fig. 1. These charac 
teristics are similar to those occurring at transonic Mach Num 


bers.'” 


DISCUSSION 


— RESULTS OBTAINED for both models are quali- 
tatively similar throughout; however, in all cases 
the flow changes introduced by wave detachment were 
found to be more violent (i.e., occur over a smaller 
range of model angles) for the two-dimensional wedge 
flow. It was found that the important parameters 
governing the position and growth of the subsonic field 
were S2/t and S;/t. It is noted that the parameter 
S./t (Fig. 2) decreases rapidly when the detachment 
angle is just moderately exceeded and then rapidly 
approaches a constant value. This confirms the theo- 
retical prediction of Busemann* and others that the 
detached flow is, in fact, relatively insensitive to the 
conditions existing at the nose of the body (i.e., nose 
angle) so long as conditions are sufficiently far removed 
from the critical angle (defined as the body angle that 
just produces sonic conditions at the body surface). 


Received February 23, 1953. 
* Now Research Engineer, Ford Scientific Laboratory, Dear- 
born, Mich. 


Hence, the important quantity governing the position 
of the subsonic field is the location of sonic velocity in 
the flow, and flow disturbances occurring near the 
position of sonic velocity (model shoulder) are thus 
relatively more important than similar disturbances 
occurring near the position of zero velocity (model 
nose). 

The parameter S,/t was found to increase through- 
out the whole range for both models used (see Fig. 2 
The continued increase at the larger angles, from what 
has been said, must then be attributed to the move- 
ment of the model tip rather than a large movement 
of the bow wave itself. It was found that the lateral 
the field, 
rapidly just as the critical angle is exceeded and then 


width of subsonic S3/t, always increases 
soon approaches a constant value. 

In conclusion, it was found that beyond the critical 
angle the subsonic field rapidly shrinks in the axial 
direction but grows, also rapidly, in the lateral direc 


tion. A stable position is approached for both models 
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somewhat above the critical angle, so that no large 
change in the location of the subsonic field occurs in 
this region. This approach is more rapid in the two- 
dimensional case. 

Along any given stream line in the subsonic field, it 
was found that the fluid particles encounter two dis- 
tinct regions in which flow gradients exist. The mag- 
nitude of these gradients is greatest at the model sur- 
face and is rapidly attenuated at moderate distances 
irom the model surface. Downstream from the bow 
wave, the fluid particles are first decelerated through a 
region of rising pressure. Next, the stream-line par- 
ticle arrives at a position where the fluid pressure at 
tains a maximum value. Finally, a region of falling 
pressure is encountered as the sonic line is approached. 

From the with the 
wedge model, it was found (as had been theoretically 
predicted by Busemann) that the detachment process 
isa smooth and continuous one except in the neighbor- 
hood of the model tip, where flow changes were found 


pressure measurements taken 


to be violent (a much reduced velocity gradient exists 
during detachment at the cone tip). In fact, from the 
measurements that were taken, it was found that, at 
all positions of the model chord except the tip, the 


flow changes introduced during wave detachment (and 


beyond) are less severe than those occurring just be- 
fore the critical angle is reached. It was found from 
the plotted distribution curves that the velocity in the 
immediate vicinity of the model leading edge, just as 
detachment begins, could be qualitatively expressed by 


Ys == K ag 


(see Fig. 1), where A, and 7 are positive constants with 
n less than unity. This type of velocity distribution is 
seen to be equivalent to that found in a sharp potential 
concave corner of angle (180 — 26s). Thus a simple 
analogy to the flow existing at the wedge tip, just at 
detachment, is possible. This qualitative formula 
indicates that the flow singularity at the tip possesses 
an infinite velocity and pressure gradient. The reason 
such a strong singularity may exist in a real fluid is 
due (it is believed) to the fact that the effect of the 
boundary layer at the extreme tip will be extremely 
small. Conceivably, the wedge angle may be altered 
slightly by the growth of the boundary layer, but the 
stream lines in the immediate vicinity of the tip must 
still retain an extremely small radius of curvature. 
Strong pressure gradients were also measured just 
ahead of the wedge corner once the flow became de 


tached. 
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The pressure distributions obtained during wave 
detachment with the cone model are shown in Fig. 3. 
The actual detachment angle for the cone model is of 
some interest, since a theoretical method of determin- 
ing this angle does not exist at present. From the 
measurements, it was found that wave detachment does 
not occur until an angle slightly greater than 6, is 
attained (64 being the maximum cone angle permitting 
an attached shock based on an infinite model chord). 
Thus detachment occurred at 46°, while 64, was calcu- 
lated as 45° 40’. 
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indicate that 
conical flow (constant surface pressure) is maintained 


The pressure distributions (Fig. 3) 


until sonic velocity is attained at the model surface 
(i.e., when the critical angle of 42° 30’ is reached 
For all angles greater than the critical, the flow ahead 
of the model shoulder is affected by the presence of the 
shoulder and the effects of a finite chord are then felt. 
The two end points of all pressure curves (Fig. 3 
were obtained on the basis of the following assumptions 
It was assumed that the pressure at the very tip would 
be identical to that on a model of infinite chord 
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ie., the expansion at the corner would not appreciably 
alter the pressure at the tip. (Since wave detachment 
js not complete throughout the range of angles shown, 
the tip does not experience stagnation conditions. ) 
At the model corner it was assumed that sonic velocity 
was attained. The ‘‘fit’’ of the end points of the curves 
calculated in this manner with the remaining experi- 
mental points is seen to be good in all cases. It is noted 
that, as soon as the critical cone angle is exceeded, the 
shoulder of the model experiences decreasing pressure 
as the cone angle is increased (because of the increasing 
bow-wave curvature). The remainder of the model 
chord, however, experiences slowly increasing pressure 
as the detach:nent angle is approached. It is seen that, 
with the growth of the subsonic field, the deviations 
from conical flow become more pronounced at all posi- 
tions of the model chord. 

The pressure measurements taken were found to be, 
on the whole, relatively unaffected by the associated 
model boundary layer. It was found, however, that 
this layer can play an important role in the viemity of 
the model corner. When no favorable potentiai pres- 
sure gradient acts on the boundary layer ahead of the 
wedge corner (i.e., before the critical angle is reached), 
the effective flow corner was moved ahead approximately 
15 boundary-layer thicknesses. . This effect is attrib- 
uted to the rounding off of the geometric corner by 
the boundary layer and to the upstream propagation 
through the boundary layer of the low pressure en- 
countered just downstream from the corner. Concern- 
ing this latter effect, little is known at the present time, 
so that even a rough check on the propagation distance 
cannot be made. It is known that moderate compres- 
sion waves can be propagated upstream through a 
laminar boundary layer for distances equal to 20 times 
(For the model used, the boundary 
Once 


the above value.* 
layer at the shoulder was certainly laminar.) 
detachment occurred, however, this difficulty was 
avoided, and a pressure close to the actual sonic value 
was measured at the This fact is 
apparently explained by the development of a favorable 
potential pressure gradient on the boundary layer at 


model shoulder. 


the wedge corner by wave detachment. 

The forebody pressure drag coefficient has been 
determined for both models in the neighborhood of the 
critical angle (Fig. 4). For both models this coefficient 
rises at an increasing rate as detachment approaches. 
At the detachment angle no sudden increase is obtained; 
rather, this angle appears as an inflection point on the 
dragcurve. Thus the high-pressure region encountered 
at the tip of both models is confined to such a small 
portion of the model chord that no discontinuous drag 
rise occurs. For the finite cone model, a considerable 
reduction in drag coefficient is attained over that pre- 
dicted for an infinite cone as soon as conical flow is 
destroyed. 

The solution used for the detached flow about a 
wedge model is carried out in the hodograph plane. 
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In this plane the boundary-value problem as originally 
proposed by Frankl’ and Guderley® has been solved 
numerically. The stream function for the detached 
subsonic flow was found as a power series of the type 


¥(r, 0) = >A, sin 2n6r"y,(7) 
n 


where 
T = (local velocity /limit velocity)” 
6 = inclination of velocity vector to wedge axis 
yn(7) = hypergeometric integral 
A, = constants determined from the flow bound 


ary conditions 


The necessary constants A, were obtained numerically 
from boundary conditions along the shock polar and 
along the supersonic expansion at the wedge shoulder. 
The supersonic boundary conditions are related to the 
subsonic-flow solution along the sonic line by means of 
an equation given originally by Tricomi.’ This solu- 
tion is only an approximate one, since only a limited 
number of constants A, are determined for a corre- 
sponding number of prescribed boundary conditions. 
However, the accuracy may always be improved if 
required. The pressure distribution at the wedge sur 
face calculated from a typical application of this method 
(using an eight-term expansion for ¥) was found to be 
in good agreement with that measured experimentally 
(Fig. 5). Further details on this entire problem may 
be obtained from reference 8. 
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The Upwash Correction for an Oscillating 


Wing in a Wind Tunnel’ 


THEODORE R. GOODMANT 
Cornell Aeronautical Laboratory, Inc. 


SUMMARY 


The upwash correction has been calculated for an oscillating 
three-dimensional wing in a circular wind tunnel at incom- 
pressible speeds. The upwash is given both at the wing itself 
wd at downstream positions where, conceivably, a tail may be 
located. The tunnel upwash is calculated numerically for the 


entire frequency range 


SYMBOLS 


p perturbation velocity potential 
frequency of oscillation 
uy rectilinear coordinates (x in direction of free stream, 


y in spanwise direction ) 


u v component of perturbation velocity in zs = 0 plane 
typical tunnel dimension 

I circulation 

t reduced frequency = wl/V 

V = free-stream velocity 

2s, 2c wing span and chord, respectively 


acceleration potential 


p free-stream density 


A wing area 

C; lift coefficient [= lift/(1/2)pV2A | 

r = radial dimension (= Y y? + sz?) 

] = eylindrical coordinate (= tan~! y 

é = downwash parameter [= (O0¢%*/0s)(2/?/AC; 


Note: Other symbols of less importance are explained at the 


point where they are introduced. 


INTRODUCTION 


-— PROBLEMS WITH WHICH THIS PAPER is concerned 
have been attacked by two other writers to the 
author’s knowledge. Jones! has solved the case of a 
center of a 
two- 


finite-aspect-ratio wing placed at the 
rectangular tunnel; MReissner® has solved the 
dimensional case. 

The present report deals with a geometrical tunnel 
configuration hitherto unexplored—namely, the tunnel 
of circular cross section. It is a lifting wing that 
experiences the interference. In addition to finding 
the upwash correction at the wing itself, the upwash 
is also calculated at various tail positions downstream 
of the wing. The problems are solved under the as- 
sumptions of incompressibility and linearization of the 
equations of motion. The author has superficially 
explored the effects of compressibility, but the results 
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are expressed in terms of untabulated functions of two 
variables. Moreover, the upwash downstream of the 
wing itself in a rectangular tunnel also involves un- 
tabulated functions of two variables even for incom- 
pressible flow. However, for the geometrical con- 
figuration with which this report is concerned, the 
upwash may readily be computed at the wing and at 
downstream positions. 

One may ask the question: Under what conditions 
is it necessary to apply wind-tunnel wall corrections to 
an oscillating wing? The author has solved elsewhere’® 
the case of the lifting wing midway between two 
vertical walls and also midway between two horizontal 
In fact, the present paper is a condensation of 
reference 3. Now if the results of these two cases, as 
well as the present case (the circular tunnel), are con- 
sufficient 


walls. 


sidered to exemplify such corrections, a 
condition for applying these corrections is obtained in 
general—namely, if the corrections would not have 
been applied in the steady case, then for the nonsteady 
case of congruent geometry they need not be applied. 
The converse is not necessarily true, however. Again, 
taking the three cases as typical, it is found that, for 
the high frequencies encountered in flutter testing, the 
corrections are negligible, but for the low frequencies 
encountered in stability testing, the corrections are of 
the same order of magnitude as the steady corrections. 

It should be emphasized that the reduced frequency, 
w// I” encountered in wall corrections in general involves 


a length, /, which is a typical tunnel dimension (e.g., 


tunnel radius) and not a wing dimension (e.g., wing 
chord). 
Among aeronautical engineers, it may be well to 


point out that the corrections here obtained involve 
the same assumptions (aside from nonsteady flow) as 
those obtained by Glauert.* 


(GENERAL CONSIDERATIONS OF ANALYSIS 


For incompressible flow, the velocity potential ¢ 
Now let ¢ be 


satisfies Laplace’s equation, Ag = 0. 
the potential of the wing in free air, and let ¢* be the 
potential of the additional flow appearing because of 
the tunnel walls. Then, Ag* = 0. There must be 
no normal flow at the tunnel wall so that 


(0/On)(o + o*) = 0 (1) 


where m is the normal to the wall. 
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The first problem at hand is that of finding an 
expression for ¢ in the case of the oscillating wing which 
corresponds to a horseshoe vortex of zero span in the 
steady case. Since this is not general knowledge, a 
brief derivation will now follow. 

Let ¢ = o(x, y, Je; then Oy,/dx may be repre- 
sented according to lifting-surface theory by [see, for 
example, Eq. (57), reference 5] 


Ov, l [Fi uoz didn (9 
= Z) 
Ox or [(x = £)? a (y = n)* + 27)" 2 \ 


where integration extends over the wing and the wake 
and all lengths are made dimensionless by dividing 
From reference 5 
In the wake 


by /, a typical tunnel dimension. 
the following equations are also obtained. 


GQ = (1 2)T (ve ates + (I (3) 
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my, 
(1/2)r, = / Uy AX (4 
x 


l 


where 


and x,, x; are equations of the trailing and leading 
edges, respectively. 

Now assume wu is constant and that the planform js 
y Then from 


rectangular (chord = 2c; span = 23s). 


Eq. (4), I, is constant and equal to 
r. = 4nu (5 
or 
Uy = lr, tc ({) 
Differentiating Eq. (3) with respect to x, the expres. 
sion for wp in the wake is 


as, Se ; . 
uy = —tkT,e°™™ ? (7 


Upon substituting Eqs. (6) and (7) into Eq. (2), there is obtained the following expression for Oy; Ox in space. 


Ov, sl", # le J dn 
= rf ce - - pe 2 
Ov S1rc , . ‘ [(x — §)? + (y — n)? + 32)’ 


tké ns dn 
ee dé ; a 
soot eS Fy er 8 


ikz ie 
re / 
tar Jt 


(S 
As s > O and c > O, this reduces to 
Ov, sT 2 j l ‘k f é - dé ’, 
= < : — 1R (9 
Ox Qe U(x? + y? + 27)” Jo [vw —-é?+yvV4+ 2 f 


From this, it is readily ascertained that 


o*y . OY —3sT',x 
es 
Ox” Ox 25x" -> y* =P 8° 
Upon integrating with respect to x from —© to 
vy, there is obtained 
Oy, ‘ sT',.2 
+ ky, = Ee oe CiT) 
Ov Qa(x* + y? + 27) 


Now define the acceleration potential, y, such that 
Oy, 


(x, y, 2) + 1Re,(x, y, 2) 
Ox 


¥(x, ¥, 3) = 

Solving this first-order differential equation with the 

condition that y, vanishes at «x = — ©, there is obtained 

“ike . ‘ 

/ e “W(x — & vy, 3) dé (12) 
0 


os ogee 


ge = 
The acceleration potential for the oscillating doublet 
is clearly 


wv = sTys/2a(x? + y? + 22)" (13) 


It will be recognized that this expression is identical 
to the steady-state acceleration potential for a horse- 
shoe vortex of zero span. If an acceleration potential 
y* corresponding to ¢* is introduced, then it is easily 


seen that Ay* = 0 and the boundary condition, Eq 


(1), is also valid for the acceleration potential, 
(0/on)(y + y*) =0 (14 


Hence, in order to solve the oscillating case, merel) 
find the acceleration potential for the steady case oi 
and then convert to velocity 


congruent geometry 


potential via Eq. (12). 

Before proceeding to particular cases, it is necessary) 
to find an expression for IT’. in terms of the lift. By 
utilizing Eqs. (81) and (S86) of reference 5, as well as 
Eq. (6) of the present report, there is found 

lift = 2p)7/*sT.(1 + tke (15 

With the conventional definition of lift coefficient, 

there is found 


ST, = AC,/4P?(1 + tke (16 


where A is the area of the wing. Asc—0, this becomes 


sT, = AC, }/? (17 


This is the well-known Kutta-Joukowski theorem for 
lift and is valid only to the first order in the chord in 
the nonsteady case. 

It should be pointed out that the assumptions that 
the plan form is rectangular and wu constant are not 


essential to the results thus far obtained. The charac 
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ter of the planform shape and up distribution are lost 
when s and ¢ are allowed to approach zero. 


THE CIRCULAR TUNNEL 
Fortunately, the steady-state acceleration potential, 
y*, has already been obtained by von Baranoff® by 
Fourier integral analysis. The solution he obtains is 
(the notation has been made to conform with this 


rep wrt) 


Substituting Eq. (18) into Eq. (12), there is found, for t 


sT’, cos 0 - f(g 
y* = = cos qxI(qr)° dq (18 
0 q 


= 
where 

fig) = (qKilqg) — g?A2(q)) Lhi'(¢ (19 
The Bessel functions appearing in Eqs. (18) and (19 


are in conformity with the definitions given by Watson.’ 
Also, the length, /, is taken as the tunnel radius. 


he velocity potential, 


sT. cos 0 “a F = f(q) 
x = , - | e™ d I cos g(x — §&)I(qr)° d dq (20 
i 0 0 q 


The radial component of velocity is 


O¢,.* sT,cos@ {[* ee i aise ; 
= = e dé cos g(x — &)I"(gr)f(q) dq (21) 
or 1* 0 0 


ee kt ™ . i . ‘ ‘ ‘ . 
Writing e '* cos RE — i sin RE and after some manipulation, there is obtained, at the center of the tunnel r = 


0, the upwash, as follows: 


ae ot, - a 7 
= i dé - cos [gx + &(k — q)\|f(q) dg + 
Oz tr- Jo l. 0 : 


| cos [gx — &k + g)|f(q) dq — 2 I sin [gx + &(R 
0 J0 


te 
i) 


— g)|f(q) dq -—1 [ sin [E(k + gq) — gx] f(q) ag (23 
0 


Now each of these integrals may be simplified according to the Appendix, and the result is 


O¢,* Oi : “= g sin gxf(q) d “© cos gxf(q) dec é 
ae k f(k)e~ i +f ae i I PIG) CF (23 
Oz 2a* 2° 0 Gr = k? 0 = k? 
Let sP, = AC,/4/? according to Eq. (17) and define (see reference 4) 
6 = (O¢,*/0z2)(l?/AC,) (24) 


Then it is found that 


1 | f(R) | (° qsin gxf(q) d ik {°° cos gxf(q) d 
- : onthe | q sin gxf(q) dq ‘4 / cos q f q) “] (93) 
Ss 2 rT JO gq? — k? T 0 gq’ — k? 


at the wing (x = Q), this becomes 


7 l Ki 4 ik a f(q) 4 (26 
8 2 r Jo gq? — k? 


For purposes of numerical integration, it is convenient t 


aaa | fbyen™ + f = gxlfia) — 
8 r Jo g? — k? 


oe l Kk ik 
~ 8 2 il wr Jo 


o rewrite Eqs. (25) and (26), respectively, as follows: 


f(k)| dq nq ik f cos axIf(q) “ {ON (27) 
r Jo gq’ — ke 
[f(q) — f(R)] a) (98) 
q° —_ k? 
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Complex upwash (86) at wing and tail positions over the frequencey spectrum. 


There is a theorem derived by Riemann’ from which ’ Watson, G. N., A Treatise on the Theory of Bessel Function 
2nd Ed.; The Macmillan Company, New York, 1945. (I 
particular, see pp. 172, 416, and 537-550.) 

‘Franklin, Philip, A Treatise on Advanced Calculus, py 
graph in which the abscissa is the real part of 6 and the — 49-481. John Wiley & Sons, Inc., New York, 1940 
ordinate the imaginary part, then, for large x, 6 would 9 Whittaker, E. T., and Watson, G. N., A Course of Moder 
appear as a circle with radius f(k)/8. Such a graph Analysis, 4th Ed., pp. 73-74; Cambridge University Press, 1940 
is given in Fig. 1 for a useful variation of x and k. 
Instead of 6, however, 86 is charted so that, for x = 0 
and k = 0, the value is unity. This value is in agree- This section i, devoted to the simplification oi 
ment with that obtained by Glauert‘ for this integrals which are of the form 


it can be shown that as x > , the integrals appearing 
in Eq. (27) vanish. Hence, if one were to make a 


APPENDIX 


case. 


7 


ef ® fsin (oe 2 REY) g 
dé (g) d 
REFERENCES I ° J Veos (q + ket I(q) dq 


1 Jones, W. Prichard, Wind Tunnel Interference Effect on the A theorem on Riemann integration® justifies the 
Values of Experimentally Determined Derivative Coefficients for writing of this as 
Oscillating Aerofoils, British R. & M. No. 1912, 1943. 
* Reissner, Eric, Wind Tunnel Corrections for the Two-Dimen- —_ en dt * §sin (q + Rel f(q) d 
P ean “ p e ; ; * ‘ 4 dé ( ad 
stonal Theory of Oscillating Airfoils, Cornell Aeronautical Labora 20 J0 7 0 lcos (q + kes J\q : 
tory Report No. SB-318-S-8, 1947. 
3Goodman, Theodore R., The Upwash Correction for an 
Oscillating Wing in a Wind Tunnel, Cornell Aeronautical Labora- the integration on £ performed. There is then ob 


The order of integration may be interchanged and 


tory Report No. AD-744-W-1, July, 1951. tained 
4Glauert, H., The Elements of Aerofoil and Airscrew Theory, 
pp. 189-198; The Macmillan Company, New York, 1944. .. f(g) dq = 7 dq 
be Reissner, Eric, On the General Theory of Thin Airfoils for fen / a Se 
Non-Uniform Motion, NACA TN No. 946, 1944. et J0 pe? + (gq + &)? 


Svon Baranoff, A., Zur Frage der Kanalkorrektur bei kom- 
pressibler Unterschallstromung, Deutsche Luftfahrtforschung, 
Forschungsbericht, Nr. 1272, 1940. (Also appears in English this becomes 
translation as: Tunnel Correction for Compressible Subsonic 
Flow, NACA TM No. 1162, 1947.) (Concluded on page 406) 


If g + k appears in the numerator, then, as u > 0, 
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Stress Analvsis of Sheet-Stringer Panels with 
Cutouts 


HARVEY G. McCOMB, JR.? 
Langley Aeronautical Laboratory, NACA 


ABSTRACT 


A method is presented for determining the stresses in axially 
loaded sheet-stringer panels having cutouts. The method is 
based on the use of so-called ‘‘unit perturbation”’ solutions that 
ire appropriately superposed to give the effects of the cutout on 
the stress distribution that would exist with no cutout The 
analysis of stringer reinforcements and doubler plates is carried 


out in a similar fashion. 


SYMBOLS 


A = cross-sectional area of stringers 

B = (E/G) (A /bdt) (b?/L?) 

E = Young’s modulus of elasticity 

G = shear modulus of elasticity 

L = distance between ribs 

P;; = stringer load at station 7j (the intersection of 


stringer j and rib 7), positive in tension 
P;(k, 1) = stringer load at station 7, 7 due to a concentrated 
unit perturbation load at station k, / 


P;;(k, = stringer load at station 7, 7 due to a unit shear per 
turbation load about shear panel k, / 

P;;\k, 1| = stringer load at station 7, 7 due toa distributed unit 
perturbation load on stringer / between rib k and 
ribk + 1 

b = distance between stringers 

i, 7 = indices for stringers, ribs, and shear panels; 7 hori- 
zontal and j vertical 

B2 = general integers 

p = concentrated perturbation load (lbs.) applied at 
station k, / 

qi = distributed perturbation lc d (lbs. per in.) applied 
to stringer / between ribs & and k + 1 

= shear perturbation load (Ibs. per in.) applied about 

shear panel k, / 

4%; = shear flow in panel i, 7 

4:;(k, 1) = shear flow in shear panel 7, 7 due to a concentrated 


unit perturbation load at station k, / 
q;;|k, /| = shear flow in shear panel 7, 7 due to a unit shear 
perturbation load about shear panel k, / 
ij. k,/{ = shear flow in shear panel 7, 7 due to a distributed 
unit perturbation load on stringer / between rib 
kandribk + 1 
thickness of sheet 


t’ = thickness of the additional portion of a reinforced 
shear panel 

P = external load per stringer applied to a tension 
panel 


Presented at the Structures Session, Twenty-First Annual 
Meeting, I.A.S., New York, January 26-29, 1953 

* The 
members of the Structures Research 


author acknowledges the assistance of the following 


NACA: Dr. 
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65 = Kronecker delta, takes the value 1 when j = / and 
takes the value 0 when 7 # / 


y(x) = [-—3B + 2(cosx — 1)]/[8B + (cos x — 1)] 

A(x) = cosh™ |7(x) 

Mi; = Lagrangian multiplier 

x* = complementary energy 

Pi; = Poa fF 4Pyy + Pea, ; 

A;;, 4j; = operators indicating the second difference in the 7 


and j directions, respectively 


INTRODUCTION 


oo CONCENTRATIONS occur in the neighborhood 
of cutouts in semimonocoque structures. Analysis 
of these stress concentrations is beyond the scope of 
elementary theories, and the development of more 
powerful methods has been in progress for a number of 
years. The solution of the cutout in a flat sheet 
stiffened longitudinally with stringers and transversely 
with ribs is of interest not only for the stress analysis 
of such simplified structures but also for the fact that 
it may indicate methods of analysis for the cutout in 
more general types of semimonocoque structures. 

An approximate solution for cutouts in axially loaded 
skin-stringer panels is presented in references | and 2. 
This solution is based on the shear lag theory of Kuhn. 
In reference 3, this type of analysis is extended to the 
case where the coaming stringers, the stringers border- 
ing the cutout, are reinforced. In reference 4, a solu- 
tion is obtained for a cutout in the center bay of a 
three-bay panel under plane loading. 

This report presents a different approach to the 
stress analysis of cutouts in flat stiffened sheet under 
axialload. The fundamental idea is that the final stress 
distribution in the structure is a superposition of (a) the 
stress distribution that would exist in the structure if 
there were no cutout and (b) the stress perturbations due 
to the cutout. This concept was suggested by Cicala.° 

The method can account for the effects of reinforcing 
the stringers bounding the cutouts—the so-called coam- 
ing stringers——and of reinforcing shear panels around the 
cutout. 


Basic ASSUMPTIONS 


The structure to be considered is illustrated in Fig. 1. 
It consists of a thin flat sheet stiffened by stringers in 
longitudinal direction and by ribs in the transverse 
direction. The cutout is located in the center of the 


panel. It removes an arbitrary number of shear panels 








388 JOURNAL OF THE AERONAUTICAL SCIENCES 





RIBS STRINGERS 








THIN SHEET: <q L k- 
L,= 





b 
| pee 
aa 
€ < 
ee 
eee 
ales 
a 
 eenaains 
A — 
ee 
a 
Fic. 1. Cutout in skin-stringer panel. 
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(c) DISTRIBUTED PERTURBATION 


Types of perturbation loads. 


NACA 


Fic. 2. 


and interrupts the corresponding number of stringers 
and is not limited to one bay in length. The loading is 
axial and uniform at the ends of the panel. The fol- 
lowing assumptions are made: 

(a) The ribs prevent transverse displacements of the 
sheet, but they offer no restraint to longitudinal dis- 
placements of the stringers. 

(b) The stringers have zero bending stiffness and 
finite axial stiffness. 

(c) The cross-sectional area of the stringers and the 
thickness of the sheet are constant. 

(d) The sheet carries only shear stresses, and the 
stringers carry only direct stresses. 

(e) The shear flow in each shear panel is constant 
throughout that panel. 

As a result of assumption (e), the variation of direct 
stress in the stringers is linear between rib stations. 
The relationship between this structural model and an 
actual panel will be discussed later. 


METHOD OF ANALYSIS 


The final stress distribution in a sheet-stringer panel 
with a cutout can be thought of as a superposition of the 
stress distribution that would exist in the panel with- 
out the cutout and certain perturbation stress distribu- 


JUNE, 1953 


tions that arise because of the cutout. 
in the analysis, then, is to consider the structure as 
though it did not have a cutout 
The stress distribution that exists in the basic struc. 
ture because of the external loads will be called the 
The second step in the 


The first Step 


the ‘‘basic’? structure. 


‘‘basic’”’ stress distribution. 
analysis is to superpose perturbation stress distribu- 
tions on the basic stress distribution in the basic struc- 
ture in such a way as to annihilate the effects of that 
portion of the basic structure which lies within the cut- 
out. 


The Unit Perturbation Solutions 


Three types of perturbation stress distributions pro- 
vide the necessary tools for this method of analyzing 
rectangular cutouts in sheet-stringer panels with or 
without reinforcement. These perturbation stress dis- 
tributions are: (1) the stress distribution caused by a 
concentrated force, p = 1, applied to a stringer at its 
intersection with a rib and acting in the direction of the 
stringer, Fig. 2(a); (2) the stress distribution caused by 
a uniformly distributed shear flow, g = 1, applied to 
the stringers bordering one shear panel in such a way 
as to cause shear in the panel, Fig. 2(b); and (3) the 
stress distribution due to one uniformly distributed 
shear flow, g = 1, applied to a stringer between two 
ribs, Fig. 2(c). For identification, these three types of 
loadings will be called, respectively, the concentrated 
perturbation loading, the shear perturbation loading, 
and the distributed perturbation loading. Note that 
the shear perturbation stress distribution is simply a 
superposition of the stress distribution due to two dis- 
tributed perturbation loads, one on either side of a shear 
panel. 

Solutions have been obtained for these perturbation 
loadings applied to a sheet-stringer panel that extends 
to infinity in both directions. Derivations of the stress 
distributions due to concentrated and distributed per- 
turbation loadings are given in the Appendix, and the 
shear perturbation solution is obtained by superposing 
two distributed perturbation solutions. The pertur- 
bation solutions depend on one parameter only, de- 
fined as B = (F/G) (A/bt) (b?/L*). The perturbation 
solutions for values of B equal to 3, 5, 7, 9, and 11 are 
given in Tables 1 to 5. Part (a) of each table contains 
the values of P;; and g,,L due to a concentrated unit 
perturbation load p = 1 at station (0, 0), where P,; 1s 
the stringer load at station (7, j)—-the intersection of 
stringer 7 and rib 7—and q;; is the shear flow in panel 
(1, 7) the shear panel having station (7, 7) in its lower 
left-hand corner [see Fig. 2(a)]. Part (b) of each table 
gives the values of P;,/L and q;; due to a distributed 


unit perturbation g = 1 between stations (0, 0) and 
(1, 0) [see Fig. 2(c)];_ part (c) of each table gives P;; L 
and q,; due to a unit shear perturbation g = 1 on shear 


panel (0, 0) [see Fig. 2(b) ]. 
The positive direction of the concentrated, dis- 
tributed, and shear perturbation loads are as shown in 
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Perturbation Distribution Due to 
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Station (0,0) Shear Panel (0,0) 
R=3 R=3 
Stringer Loads: 4 Stringer loads: “4 
os i 
rer : = - 7 i oe 1 2 ae. L 
c <= n.16e nna 0.97 1 -7.3162 1188 ~0.0L80 ~0.9212 
1 oi ae i Se eae corte | ech | sesg | 3:0835 
0.9216 9.933 0.0390 0.0393 3 0.0063 -0.0145 “0.0163 “0.01L6 
0 9.9101 0.0181 0.0232 9259 h 4.0022 1.005 0.6079 6.008h 
058 0.0109 0.0149 0.0177 S 7.0619 0.0029 -0.00L2 0 ..COL9 
0.0038 0.0073 9.01 3 126 & —O,.cNNG ~) C016 ~0.002) 1.0030 
0.0027 0.0052 0.0075 0.9L 7 1.0903 0.0910 0.0015 —.0020 
0.0020 0.0939 0.0957 1.0073 a 0.0992 1.0006 =9,9910 0.9913 
0 0.16 9031 0.00L5 9.0059 9 91 6 .000L 0.0997 =~ ,.0009 
rn 0.9013 9.125 0.0037 ALAR 10 0.09%] 9.0903 0.0005 “0.0907 
0 9,001] 0 0021 9.0031 C eee ll “7.91 -7,.0902 -0.0ML 0.0005 
fal 0.9909 0. 1.00027 now denials pies peating - x 
0 0.0008 iors 0.93 | on91 +2 0.0002 “0.0993 -0.000L 
Panel Shears 44 4 
Panel Shears: ‘ij 
7 i=o I ae a 
— oil i 
| 9 0.12h6 0.9521 0.92L3 12° a 
| 1 9.0672 WR 9.0317 0.0221 4 1 2 3 
? 0.0356 2.031 | 0.0270 | 0.9719 0 28h 0.16ER nine | a.naeé 
i 0.0197 0.9190 0.0178 0.0164 1 =).0.933 0316 9.0017 NANAK 
Ps 0.0159 0.9156 0.919 0.01h0 2 1.9230 1.0173 0,908 —,.W22 
4 0.0132 0.0130 0.01764 9.0120 3 9.9103 “,9092 0.0056 = .9039 
7 0.9112 9.9111 0.010 g h —D.0089 0955 -0.00L6 +.003L 
2 0096 0.0096 n, S 0.0038 -).9036 0.0032 0.0027 
| 9 0.098h 0.83 0. 6 0.0027 ~).0026 0.002), 0.0021 
} 10 0.0073 0.0073 0. 7 ~).0920 —),0020 0.0919 ~.0017 
1 0.006L 0.0064 0. R .016 ~),0016 0.0015 0.091 
12 0.0056 0.0056 a 9 4.9013 ~.0013 ~.0012 -0.0012 
10 “0.0011 0.0010 0.0010 0.0010 
11 0.0099 0.0009 0.9904 ~0.0008 
TABLE 1(b) 1? —O.0R =—1,.000R8 ~0.90048 ~0.0908 
Perturbation Distribution Due to TABLE 2(a) 
Distributed Unit Perturbation Load On 
Stringer O In Bay 0O Perturbation Distribution Due to 
Concentrated Unit Perturbation Load at 
B = 3 Station (0,0) 
B=5 
Stringer Loads: Pi3 
I Stringer Loads: Pij 
3 isl 2 3 L 
r i 0 1 2 3 L 
0 0.3578 0.1920 0.1195 0.0839 Jj 
1 0.016 0.0732 0.971L 0.0627 
2 0.0115 0.0288 0.0372 0.039L ° 0.5000 0.2913 0.1853 0.1283 0.0955 
3 9.0062 0.01L3 0.0209 0.0248 1 0 0.0590 0.0738 0.072h 0.0662 
3 0.005 -C . VU. 2 0 0.0175 0.0297 0.036 0.0390 
L 1.0029 0.008h 0.0130 0.016 ° “ 
Oe - os 3 0 0.0080 0.01h8 0.0200 0.0236 
5 0.0019 0.0056 0.9086 060125 u 0 0.005 0.0087 0.0123 0.0152 
6 0.0013 0.000 0.006 0.0085 5 0 0.0030 0.0058 0.0083 0.0105 
7 0.0010 0.0030 0.90L8 0.0065 6 0 0.0021 0.00h1 0.0060 0.0077 
8 9.0908 0.0023 0.0038 0.0052 7 0 0.0016 0.0031 0.00L5 0.0059 
9 9.0006 0.9019 0.0031 0.0013 8 i) 0.0012 0.002 0.0036 0.0047 
10 0.0905 0.0916 0.0026 0.0936 9 O 0.0010 0.0020 0.0029 0.0038 
11 0.090 0.001h 0.0022 0.0031 7 0 0.0008 0.0016 0,002h 0.0032 
[12 0.0004 0.0012 0.0019 0.0027 - 0 0.0007 0.001L 0.0021 0.0028 
1 fe) 0.0006 0.0012 0.0018 0.002h 
Panel Shears: 44 
Panel Shears: %4 5! 
j i=0 1 2 3 0 1 2 3 
R) 
0 0.1422 0.0829 0.0363 0.9178 
1 0.0589 0.051 0.9380 0.0266 fe) 7 0.10bL 0.0530 0.0285 0.0164 
2 0.9359 0.031 0.9296 0.02hh 1 0.0453 0.0383 0.0298 0.0226 
3 0.0256 0.02h9 0.0230 0.0205 2 0.0278 0.0260 0.0232 0.0200 
4 0.9197 0.019L 0.0185 0.0171 3 0.0198 0.0192 0.0180 0.0165 
S 0.0159 0.0158 0.9162 0.01bb 4 0.0153 0.0150 0.01bb 0.0136 
a ¥ 7 . 5 0.012k 0.0122 0.0118 0.0114 
6 0.0132 0.0131 0.0128 0.0123 ? 
° ° ° ° 6 0.0103 0.0101 0.0100 0.0097 
7 0.0112 0.0112 0.0110 0.0106 7 0.0087 0.0086 0.0085 0.0083 
8 0.0097 0.0096 0.9095 0.0092 a 0.0075 0.007h 0.007 0.0072 
9 0.008 0.008 0.0083 0.0080 9 0.0065 0.0064 0.006k 0.0063 
10 0.0073 0.9073 0.0072 0.0071 10 0.0057 0.0056 0.0056 0.0055 
ll 0.006L 0.0064 0.006L 0.0062 11 0.0050 0.009 0.00L9 0.00L9 
| 12 0.0056 0.0056 0.0056 0.0055 12 0.00bh 0-00h3 0.003 0.00h3 
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TABLE 2(b) 
Perturbation Distribution Due to 
Distributed Unit Perturbation Load on 
Stringer O In Bay 0O 


B= 5 


Stringer Loads; Pij 
L 





TABLE 2(c) 
Perturbation Distribution Due to 
Unit Shear Perturbation Load About 
Shear Panel (0,C) 


B=5 


Stringer Loads: Piy 

























































































" L 
j 1 2 3 h 
) 0.368h0 0.2328 0.15h1 0.1105 4 1 2 3 h 
1 0.03L9 0.068) 0.0738 0.0695 j 
2 0.0091 0.02h2 0.0335 0.0380 
3 0.00L0 0.0115 0.0176 0.0219 1 6.3191 -~0.16Lh ~).0°03 -0.0h10 
, payer 0.0067 0.0106 0.0138 2 0.0258 | -0.0Lb3 0.003. | 0.0316 
3 0.0025 0400bl, 0.0071 0.009L 3 0.0050 | -0.0126 ~0.01°9 | -0.0160 
6 0.0010 0.0031 0.0051 0.0069 is 0.0017 -0.00L8 0.0070 “0 6081 
r 0.0008 0.0023 0.0038 0.0052 5 0.6008 | 0.0023 0.0035 | -0.colh 
8 0.0006 0.0018 0.0030 0.00L1 6 -0.CcO0oL 0.0013 —C.CC20 6 .0026 
9 0.0005 0.0015 0.002 0.003 7 ~0.0003 -0.6008 6.0012 0.0016 
” 0.0001 0.0012 0.0020 0.0028 E -0.0002 0.6005 ~0.0008 0.0011 
Ls. bya — 0.0017 0.002h 9 -0.0001 | ~0.cooh ~0.6006 | 0.0007 
= 0.0009 0.9015 0.0021 10 -0.0C01l -0.0002 ~0.000b 1.0006 
11 -0.0001 0.0002 0.0003 0.000) 
12 -0,.0000 -0,.0001 —.0002 5.0003 
Panel Shears: 44 
Panel Shears: 444 
4 0 1 2 3 
: i 0 3 2 3 
fe) 0.1160 0.0756 0.039 0.0218 j 
1 0.0L61 0.0421 0.030 0.0261 
2 0.0280 0.0270 0.02h7 0.0216 0 0.232C 0.1£12 0.6787 0.0436 
3 0.0199 0.0195 0.0187 0.0172 if 0.0698 0.0335 ~0.005h 0.0CL3 
L 0.0153 0.0152 0.0148 0.01L0 2 0.6182 -0.01°1 ~0.0093 ~0..COLS 
5 0.0123 0.0123 0.0120 0.0116 3 0.0081 ~0.0075 ~0.0060 -0.C0LL 
6 0.0102 0.0102 0.0101 0.0098 L -0.COL6 ~0.0chhL 6.0039 ~0.0032 
7 0.09087 0.0087 0.0086 0.008h 5 ~0.0030 0.0029 ~0.C0027 -0,.002h 
8 0.007h 0.0075 0.007h 0.0073 6 —©.0021 —C.6020 -0.0020 -0.0018 
9 0.0064, 0.0065 0.006h 0.0063 7 =0,0016 -0,0015 -0.601¢ -0,.001L 
10 0.0056 0.0057 0.0056 0.0056 8 0.0012 0.0012 -0.06012 0.0011 
11 0.009 0.0050 0.00L9 C.00L9 9 6.0010 0.0010 -0.0010 0.0009 
12 0.00L3 0.00LL 0.003 0.00L3 10 -0.C008 0.0008 -0.0008 -0.0008 
11 0.0007 -6.C007 -C.6007 6.6007 
12 —5.0006 -0.C006 —6.0006 sane 























Fig. 2; stringer loads are assumed positive in tension, 
and the sign convention for shear flows is the same as 
that for shear perturbation loadings. The perturba- 
tion solutions for arbitrary locations of the perturba- 
tion loadings are readily obtained from the tables by 
means of changes of indices. 


Analytical Procedure 


Unreinforced Cutout.—Suppose the basic stress dis- 


tribution is known. In a tension panel under uniform 
load, this basic distribution simply states that the 
stringer load everywhere is equal to the applied stringer 
load and that the shear flow everywhere is zero. The 
stress distribution in the panel with a cutout is deter- 
mined in the following manner: Concentrated per- 
turbation loads are imposed on the basic structure at 
stations where stringers will be interrupted by the cut- 
out, and distributed perturbation loads are imposed 
along stringers that will bound the cutout. The mag- 
nitudes of these perturbation loads are adjusted so that 
simultaneous removal of the cutout material and the 


perturbation loads would leave the structure outside 





of the cutout unaffected. Then, with this condition 
satisfied, the stresses in the structure with the cutout 
subjected to the basic loading only are precisely those in 
the basic structure due to the basic loading, together 
with the perturbation loadings. 

Consider the case where the cutout is only one bay 
long, as in Fig. 1. At each pair of stations in the basic 
structure where a stringer would be interrupted by the 
cutout, there is placed a pair of concentrated perturba- 
tions equal in magnitude and acting in opposite direc- 
Fig. 3 shows these perturbations applied to a 
It is to be 


tions. 
basic structure for the cutout of Fig. 1. 

noted that, under these perturbation loads, as well as 
under the basic loading, no shear flow is introduced into 
the shear panels of bay 0. 

At each point where a concentrated perturbation 
load is applied there is a discontinuity of stringer load 
equal in magnitude to the perturbation load. Thus, 
in the stringers that will be interrupted by the cutout, 
a distinction must be made between the stringer load 
just inside of the proposed cutout and that just outside 


of the proposed cutout. Consider the stringer load 
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TABLE 3(a) TABLE 3(b) 
Perturbation Distribution Tue to Perturbation Distribution Due to 
Concentrated tmit Perturbation Load at Distributed Unit Perturbation Load on 
Station (0,0) Stringer © in Bay O 
"eT B=7 
Stringer Loads: Pij Stringer Loads: "tu 
i 0 1 2 3 L i 1 2 3 L 
j 
0 0.£000 0.2168 0.2221 0.151b 0.121b2 0 0.2991 0.2593 0.1792 0.1314 
1 0 0.C53k 0.0713 0.0738 0.670? 1 0.0308 0.0642 0.0733 0.0722 
2 ie} €.C1S51 0.6266 0.0338 0.0377 2 0.0078 0.0212 0.0306 0.0360 
3 0 0.C068 0.0129 0.0178 0.0216 3 0.003L 0.0099 0.0155 0.0198 
L 0 0.0039 0.075 0.01¢8 0.0135 u 0.0019 0.057 0.0092 0.0122 
5 0 0.025 0.COL9 0.9072 0.0092 5 0.0013 0.0037 0.0060 0.0C82 
6 0 0.0018 0.0035 ¢.0051 9.067 6 0.0009 0.0026 0.00L3 0.0059 
7 0 0.0013 9.0026 9.0039 Arnel 7 0.0006 0.0029 0.0023 0.00L5 
8 ° C.0n1O 0.0020 9.0030 A.caho 8 0.0005 0.0018 0.0026 0.0035 
Be = pope hid 0.0017 0.026 9.0033 9 0.000L 0.0€13 0.0021 0.0029 
= sgesiv 0.c01, 9.71 0027 10 0.000L 0.0010 9.0017 0.002 
rT) ‘ replies ©.001? O.nn1R 9.0023 11 0.0003 0.0009 0.0c15 0.0021 
‘ A.cele o.0n1s N.cN20 12 0.0002 0.0008 050013 0.0018 
Panel Shears: 445 Panel Shears: 4 
i 0 4 0 1 2 3 
; 1 2 3 j 
- anne _ 0 0.3009 0.0699 9.0L00 0.0739 
: gptabee 0 0519 9.C304 O.01R6 1 0.0392 9.0366 0.0309 0.0250 
2 ae ap O.cere | 0.ceeP ? 0.0237 0.0231 0.0216 0.0195 
3 “O16! Y 066706 OomIFl 3 0.0169 0.0166 0.0160 0.0182 
j L 0.0130 0.0128 0.0126 0.0122 
c Ss 0.0105 0.010L 0.6102 0.6100 
é 6 0.0087 0.0086 0.0086 0.008L 
7 7 0.007h 0.0073 0.0073 0.0072 
i 8 0.006), 0.0063 0.0063 0.0062 
9 0.0056 0.005L 0.0C5L 0.005L 
10 0.00L8 0.00L8 0.0ch8 0.00h7 
11 0.00L2 0,002 0.00L2 0.00L2 
12 0.0038 0.0036 0.0037 0.0036 
TABLE 3(c) 
Perturbation Distribution Due to 
Ynit Shear Perturbation load about 
Shear Panel (9,() 
R«?7 
Stringer Loads: Pij 
L 
i 1 2 3 j b 
j 
1 0.3683 0.1951 0.1059 0.0592 
2 0.0230 ~0.0430 ~0.0127 ~0.0362 
3 ~0.COLL 6.0113 0.0151 -0.C162 
=| 0 | i 0.0015 ~0.00L?2 0.0063 0.0076 
7 5 0.0007 0.0020 ~0.0031 ~0.00L0 
I 0 I 2 6 ~0.000k 0.0011 0.0018 0.0023 
. : vy 7 0.0002 0.0006 0.0010 ©.colL 
—_ MASA 0 0.0001 ~0.000L 0.0007 0.0009 
Fic. 3. Perturbation loads for the cutout of Fig. 1. 9 0.0001 “0.0003 ~©.0005 0.0006 
10 ~0.C001 ~.c002 000k -0,0008 
1) 6.0001 0.0002 <.0C002 0.0093 
r2 0.0000 0.0001 0.0002 0.0003 
just outside of the proposed cutout in a stringer that 
will be interrupted by the cutout. When the magni- ae See 
tudes of the concentrated perturbation loads are ad- ‘ ear 7 : ; 
Justed to give, together with the basic loading, zero : 
<a ahi i ee ln . : " i a ) 0.2017 0.1398 0.6801 0.078 
stringer loads at these points, the portion of the struc au om Sua pap 
ee tenn wt a > ease — ; 2 0.0155 ~0.013k ~0.009L: ~).0055 
ee lying within the cutout can be removed from the 3 ae Py = <= 
ae . aes 3 Pee a — wali 0.0039 ~0.0038 ~0.003h 0.0030 
asic structure, and, simultaneously, the perturbation ; Fray Pay ae en 
loads themselves can be removed from the basic struc- 6 0.0018 ~0.0017 “0.0017 0.0016 
; , , 7 0.0013 0.0013 ~.0013 ~0.0012 
ture without affecting the rest of the structure. Thus, 8 “0.0010 “0.0010 0.0010 0.0010 
: : s : . 9 0.0008 0.0008 0.0008 ~0.0008 
in this problem of a one-bay-long cutout, distributed 10 ~0.0007 0.0007 ~0.0007 “0.0007 
i 5 11 0.0006 ~0.0006 ~0.0006 “0.0006 
perturbation loads along the coaming stringer are not 12 “0.0005 0.0005 0.0005 “0.0005 


needed. 
The conditions that the concentrated perturbation 
loads must satisfy are expressed mathematically in the 


following manner: Let /,, represent the magnitude of 


































































































Let Pi;(k, 1) represent 
the stringer load at 7, 7 toward the outside of the cut- 


the perturbation load at k, /. 


out due to a unit concentrated perturbation load at 








392 JOURNAL OF THE 


AERONAUTICAL 


SCIENCES JUNE, 1953 





f= 
= 























2 7 - 
| a” 
0 
0 & —EEE 
mn 
“5 
“2 








- sie ee . | 
— 0 | 2 

















Nata 









































Fic. 4. Perturbation loads for a cutout extending two bays in 
length. 
TABLE h(a) 
Perturbation Distribution Due to 
Concentrated Unit Perturbation Load at 
Station (0,0) 
B-9 
Stringer Loads: Pij 
i 0 1 2 3 L 
J 
0 0.£000 0.3328 0.2322 0.1699 0.1300 
1 te) 0.0493 0.0687 0.0738 0.0722 
> 0 0.0135 0.62b3 0.0317 0.6362 
3 0 0.0060 0.0115 0.0162 6.0200 
L (a) 0.03L 0.0067 0.0097 9.0123 
S 0 0.0022 0.00LL 0.6L 0.0083 
é Le] 0.0016 0.6031 0.006 0.0060 
7 0 n,0012 0.0023 0.003L o.cohe 
p 0 0.0009 0.c018 n.0027 Maas 3 
9 0 n.0nn7 0,015 9.902? 0.0029 
1° ¢ 9.0006 0.0012 9.0018 n,002h 
1 ( 0.0005 0.0010 0.0016 o.n02 
12 C 9.0005 0.0009 0.cO1L 0.0018 
Panel Shears: 9; 3! 
i © 1 2 3 
3 
0 0.0836 0.0503 0.0312 0.0199 
1 0.033 0.0309 0.0261 0.0215 
2 0.0209 0.C201 0.0187 0.C170 
3 0.C1L8 0.C1L6 0.01L0 0.0133 
b 0.011 0.0113 0.0110 0.C107 
S 0.0092 0.0091 0.0090 0.0088 
6 0.0077 0.6076 0.0075 0.cO7kL 
7 0.0065 0.0065 0.006L 0.0063 
8 0.0056 0.6056 0.0055 0.0055 
9 0.00L8 0.col8 0.CcOLS 0.cOL8 
10 0.00L2 0.COh2 0.COL2 0.00h2 
ll 0.9937 0.0037 0.0036 0.37 
1? 0.0033 9.0033 0.0032 0,32 
—_ 























k, |. Let P represent the basic stringer load. The 
conditions on the p,; are expressed as 


> Dd pePil(k, ) + P =0.... (1 
k l 


where the summation is carried over the k, / for all per- 
turbation loads and where the equations are written 
for each station 1, 7 at which the cutout will interrupt a 
stringer. 


TABLE h(b) 
Perturbation Distribution Due to 


Distributed Unit Perturbation Load on 
Stringer O in Bay O 


Stringer Loads: Pag 



































I 
i 1 2? 3 h 
a ae 
0 0.L09L 0.27AS 0.1987 0.1L 85 
1 0.020 0.0606 9.9720 0.0733 
2 0.0069 0.0192 0,0283 n 03h2 
3 0.0030 0.0089 0.01L0 0.0182 
L 0.0017 0.9050 0.0082 0.0110 
S 0.0011 0.0033 0.C0SL 0.0073 
A 0.0008 0.9023 0.0038 0.0053 
? 0.0006 0.0017 0.0029 0.00L0 
R 0.0005 9.001L 0.0023 0.0031 
9 0.000L 0.0011 0.0018 0.0025 
10 0.0003 0.0909 0.0015 0.0071 
11 0.6003 0.0008 0.0013 0.0018 
12 0.0002 0.0007 9.0011 0.0016 
Panel Shears: Wij 

a 

i 0 1 2 3 | 

3 | 

| 

0 9.0906 0.065L 0.0399 0.0251 | 

a. 0.03L7 0.0328 0.0285 0.0238 

2 0.0209 0.0205 9.019b 0.0179 | 

3 0.01L8 0.0147 0.013 0.0137 | 
h 0.011h 0.0113 0.0112 0.0109 

5 0.0092 0.0092 0.0091 0.0089 | 
6 0.177 0.0076 0.0076 0.0075 
7 0.0065 0.0065 9.0065 0.006L 
8 0.0056 0.0056 0.0055 0.0055 
7 0.00L9 0.COL8 0.00LE 0.00L8 
10 0.00L2 0.00L2 0.002 0.00L2 
ll 0.0037 0.0037 0.0037 0.0037 
12 0.0032 0.0033 0.0032 0.0032 




















In case the cutout does extend more than one bay in 
length, distributed perturbation loads must be im- 
posed on the basic structure on the stringers bounding 
the edges of the cutout (the coaming stringers) within 
the bays where the cutout exists. Fig. 4 illustrates the 
application of perturbation loads when the cutout 
extends two bays. The magnitudes of the distributed 
perturbation loads must be adjusted so that, in each of 
the shear panels adjacent to these perturbation loads 
on the inside of the proposed cutout, the sum of the 
shear flows due to all the perturbation loadings just 
cancels the distributed perturbation loading at that 
shear panel. The conditions on the stringers at stations 
where they are interrupted by the cutout remain the 
same as before. When these conditions on the dis- 
tributed and concentrated perturbation loads _ hold, 
the boundary conditions at the cutout are satisfied, and 
the portion of the basic structure within the proposed 
cutout can be removed along with the perturbation 
loads without affecting the rest of the structure. 

The mathematical formulation of the problem now 
Let P;; |k, /| mean the stringer 
load at 7, j due to a distributed unit perturbation load 
Let q,;(R, !) 


becomes as follows: 


on stringer / between rib k and rib k + 1. 
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TABLE U(c) 


Perturbation Distribution Pue to 
nit Shear Perturbation Load about 
Shear Panel (0,9) 












































R=9 
Stringer Loads: od 
i 1 2 3 u 
4 
1 —0.3Flk 06,2179 0.1267 0782 
2 0.0211 -C.0h1L 7.0137 0.0367 
3 0.0039 0.0103 0.01L3 -0.0160 
k -6.0013 6.0038 0.0058 “0.0972 
5 0.0006 -#.0018 0.0028 1.0736 
e 0.0003 “0.0010 -0.0015 0021 
7 6.0002 “0.0006 —©.0010 ~.0913 
’ 0.0001 ~0.000k 0.0006 0.0008 
0.0001 0.0002 ¢.000h —.0006 
0.0001 0.0002 0.0003 =—0.000L 
11 6.0000 ¢.0001 0.0002 —.0003 
12 —,0000 0.0001 0902 0.0002 
Panel Shears: 444 
a 
i 0 1 2 3 
J 
0 0.1813 0.1309 0.0798 0.0502 
1 “0559 0.09326 ~.01lh 0.0013 
2 0.0138 0.0123 “0.91 0.9059 
3 1.60 =0.0958 “0.9951 0.09? 
L 0.003L -0,.093L 70.0931 ~.0028 
5S 6.0072 0.0922 0.0021 =0.9020 
6 0.016 0.0015 0.0015 0.1L 
7 .0C12 7.011 0.11 ~0.M11 
Q 0.909 0.009 0.0009 0.0009 
9 0.0/8 -0.0007 0.9007 ~0,0907 
10 0.9006 0.9006 0.0006 0.0906 
il 0.0095 0.0005 0.9005 0.0005 
12 6.000L 0.0095 0.090 0.0005 























represent the shear flow in panel 7, 7 due to a concen- 
trated unit perturbation load at station k, /, and let 
gy{k, /| represent the shear flow in shear panel 7, 7 due 
to a distributed unit perturbation load on stringer / be- 
tween rib k and ribk + 1. If G, represents the mag- 
nitude of the distributed perturbation load on stringer 
! between rib k and rib k + 1, the equations that 
must hold are 


OY pu Polk. D+ LY Ges Po (b+ 
a | k l 


DD Per ule, ) + UD Gequik, = Gy (3) 
ae te ns 
where the summations over k and / are extended over 
all the concentrated perturbation loads and all the dis- 
tributed perturbation loads. Eqs. (2) are written for 
every station 7, 7 where the cutout will interrupt a 
stringer. Eqs. (3) are written for each shear panel , 7 
adjacent to a coaming stringer and inside of the pro- 
posed cutout. 

Solution of the system in Eq. (1) for a cutout one 
bay long or solution of Eqs. (2) and (3) for cutouts ex- 
tending more than one bay in length yields the required 
magnitudes of the concentrated and distributed per- 
turbation loads f,, and g,. The quantities P;,(k, /) 
and q,,(k, /) in Eqs. (1) to (3) are obtained from parts 


TABLE 5S(a) 


Perturbation Distribution Due to 
Concentrated Unit Perturbation Load at 
Station (0,9) 


B = ll 


Stringer Loads: P44 





























i 0 1 2 3 L 
J 
0 0.5000 0.3456 0,281 0.1852 9.1435 
1 0 0.0460 0.0663 0.0731 0.0732 
2 0 0.0123 9.0225 0.0299 0.0308 
3 ra 9.9055 0.0106 0.0150 9.9196 
L 0 9.0931 0.61 0.0088 9.9113 
4 0 9.9920 a.m 0.0958 9.9076 
6 0 0.001h 0.9028 0.00h2 n.05h 
7 0 9.0910 9.0021 0.0031 0.Mbh1 
8 0 0.0008 9.9016 0.902h 0.0932 
9 0 0.0007 0.9013 0.020 9.126 
10 0 0.0906 0.0011 0.17 9.0022 
11 0 0.0005 0.09 0.091h ¢.0019 
12 0 0.000L 0.0008 0.0012 0.9016 
eaamiemenanentil 





Panel Shears: % 41 























i 0 1 2 3 

j 

0 0.0772 0.0487 0.0315 0.0209 
1 0.0312 0.0285 0.027 0.0208 
2 0.0189 0.0183 0.0172 0.0159 
3 0.013L 0.0132 0.0128 0.9123 
u 0.9103 0.0102 0.0109 0.0098 
S 0.0083 0.0083 0,0082 0.0080 
6 0.0069 0.0969 0.0068 0.0067 
7 0.0059 0.0058 0.0058 0.0058 
8 0.0050 0.9050 0.50 0.0050 
9 0.MbLk 0.00bLL 0.00bL 0.90L3 
10 0.0038 0.9038 0.0038 0.9038 
11 0.003h 0.9033 0.0033 0.0033 
12 0.0029 0.0029 0.9029 0.0029 








(a) of Tables 1-5; Pi; {k, /} and qj; {k, /} are found in 
parts (b) of these tables. The stresses about the cut- 
out are then obtained by superposing the perturbation 
stress distribution in the basic structure on the basic 
stress distribution. 


The unit perturbation solutions presented in the 
tables apply to panels that extend to infinity in both 
directions; they can be applied directly to a finite 
panel without serious error if the cutout does not come 
too close to the edge of the panel. A refinement in the 
analysis to account for edge effect will be discussed later. 

To illustrate the method of analysis, the solution 
for a cutout in a tension panel will be presented. The 
case to be analyzed is that of a long panel with 16 string 
ers having a cutout one bay long which interrupts four 
(See Fig. 1.) 


Other pertinent properties of the panel are taken as 


stringers and removes five shear panels. 


follows: 


b = 3.125 in. 

L = 3.00 in. 

t = ().0331 in. 

A = 0.2516 sq.in. 


A/bt = 0.2516/0.1035 = 2.43 
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TABLE 5(b) 
Perturbation Distribution Due to 
Distributed Unit Perturbation Load on 
Stringer O in Bay 0 


B=1l 


Stringer Loads: “i 





























1 2 3 l 
j 
0 0.4169 0.2933 0.2145 0.1630 
1 0.0258 0.0576 0.0705 0.0735 
2 0.0063 0.0176 0.026L 0.0325 
3 0.0027 0.0080 0.0128 0.0169 
h 0.9016 0.0046 0.0075 0.0101 
. 0.0010 0.0030 0.00L9 0.0067 
6 0.0007 0.0021 0.0035 0.9048 
7 0.0005 0.0016 0.0026 0.0036 
8 0.000h 0.0012 0.0020 0.0028 
9 0.0903 0.0010 0.0016 0.0023 
10 0.0003 0.0008 0.0011 0.0019 
11 0.9002 0.0007 0.9012 0.0016 
12 0.0002 0.9006 0.0010 0.0014 





Panel Shears: jj 


























| 0 1 2 3 

0 0.0831 0.0618 0.039h 0.0258 
1 0.0315 0.0300 0.0266 0.0227 
2 0.0189 0.0186 0.0178 0.0166 
3 0.013h 0.0133 0.0130 0.0126 
4 0.9103 0.0103 0.0101 0.0100 
5 0.008) 0.0083 0.0082 0.0081 
6 0.0070 0.0069 0.0069 0.0068 
7 0.0059 0.0058 0.0058 0.0058 
8 0.0051 0.0050 0.0050 0.0050 
9 0.00L5 0.00hLL 0.0043 0.00Lh 
10 0.0039 0.0038 0.0038 0.0039 
11 0.0035 0.0033 0.0033 0.0034 
12 0.0031 0.0029 0.0029 0.0030 





b?/L? = (3.125/3.00)? = 1.085 
B = (E/G) (A/bt) (b?/L?) = 6.99 


beP wl, 2) + puPwe(l, 1) + poPe(l, 0) + pr, -1 Pel, 
From the conditions of symmetry, 

Px = —pu = —po -1 = pi -1 

Pu = — po ans — Poo = Pio 
With the assumption B = 7 it is found from Table 


3(a) that 


P,(1, 2) = 0.5000 

P,(1, 1) = Pe(l1, 0) = P»(1, —1) = 0 
P,,(0, 2) = 0.3158 
P,(0,1) = 0.0534 


TABLE 5(c) 


Perturbation Distribution Due to 
Unit Shear Perturbation Load about 
Shear Panel (0,0) 


B=1ll 


Stringer Loads: “i 









































‘ 1 2 3 L 
J 
1 0.3910 0.2356 “0.040 0.0895 
2 0.0196 0.9400 ~0.OhL1 ~0.9L10 
3 0.0035 0.9096 0.0136 0.0156 
L ~0.0012 0.9035 0.905L 0.9068 
5 0.0006 0.0016 0.0026 0.003 
6 0.0003 0.9009 —0.901k 0.9019 
7 0.9002 0.9005 =0.9008 0.9012 | 
8 0.9001 -0.900L 0.0006 0.9008 
9 0.0001 0.0002 ~0.900L 0.0005 
10 0.9001 0.002 0.0003 =0,.900L 
11 0.9000 -0.9001 0.9002 0.0003 
12 0.0000 -0.9001 “0.9002 -0.9002 
Panel Shears: qi 
i 0 1 2 3 
Jj 
(a) 0.1663 0.1236 0.9788 0.0515 
1 ~0.0517 0.0318 0.9128 0.0030 
2 0.0125 0.0114 0.0088 “0.0061 | 
3 0.0955 0.0053 -0.00L8 ~.00h1 
iG 0.0031 0.0030 0.0029 0.026 
5 0.0020 ~0.0020 ~0.0019 0.0018 
6 ~0.001L 0.001 -0.091h -0.0013 
7 6.0011 =0.9010 0.9010 -0.0010 
8 0.9008 0.0908 0.0008 0.0008 
9 0.9007 0.0007 0.0007 0.0007 
10 =0.0005 0.0006 0.0006 = ..0005 
11 20005 “0.0005 0.0005 0.0005 
12 0.000 ~.900L 0.900; 0.900 




















Fig. 3 illustrates the perturbation loads to be applied 
to the basic structure; the perturbations are doubly 
symmetric about the cutout. 

Application of Eq. (1) in order to liquidate the load 
in stringer 2 just outside of the cutout at (1, 2) (see 


Fig. 3) gives 


vied 1) + PoP (0, 2) +- buP (0, 1) + 


boP (0, 0) + po, 1P2(0, —1) + P = 0 (4) 
P»(0, 0) = 0.0151 
P,.(0, —1) = 0.0068 


Substituting these values into Eq. (4) yields 
P — O0.1774py + 0.0685p = 0 (5 


Similarly, liquidation of the load in stringer | at (1, 1 


requires that 


P + 0.0685py — 0.1308h,, = 0 (6 
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TABLE A 
_——Stringer Loads- -Shear Flows— 
Intersection of Rib 

and Stringer Load Shear Panel Shear Flow 
(1, 3) 1.906P (1, 0) OP/L 
(1, 4) 1.331P (1, 1) —0.176P/L 
(1, 5) 1.176P (1, 2) —0.515P/L 
(1, 6 1.1138P (1, 3) —0.074P/L 
(1, 4) —0.021P/L 


Because of symmetry, Eqs. (5) and (6) are the only 
independent equations. That is, if equations are 
written stating the liquidation of the load at any of the 
other points where a stringer is cut, they will turn out 
to be identical to either Eq. (5) or Eq. (6) 

The solution to the system of Eqs. (5) and (6) yields 


Pw = 10.75P, bu = 13.28P 


Stringer loads and shear flows in the neighborhood of 
the cutout are obtained with this solution in conjunction 
with Table 3(a). Some results are shown in Table A 
Reinforced Cutout. 
to include the problem of shear reinforcement about the 
Assume that some of the shear panels in the 


The analysis is easily extended 


cutout. 
neighborhood of a cutout are reinforced by the addition 
of a certain thickness of sheet (i.e., a doubler plate). 
Then, the procedure consists of adding shear perturba- 
tions to each of the panels of the basic structure which 
is reinforced in the actual structure. On the doubler 
plate, or the added portion of skin which is the shear 
reinforcement of the panel, place the same perturbation 
with the opposite sign. Then, for each reinforced panel, 
write an equation that states the requirement that the 
shear stress in the shear panel of the basic structure shall 
equal the shear stress in the excess portion of skin used 
to reinforce that panel. When this condition is satis- 
fied, the loaded excess portions of sheet (the doubler 
plates) can be considered inserted in the basic structure 
without disturbing continuity. The shear flows on the 
upper and lower edges of the excess portions of sheet 
cancel the shear perturbation loads on the basic struc 
ture, and the shear flows on the vertical edges of the 
excess portions of sheet go directly into the rigid ribs. 


Assuming B = 
P — 0.5000Py + 0.3158fhn + 0.0534py, 


Condition (b) is automatically satisfied by symmetry. 


- 0.3683¢nL —0.1951gnl — 0.0230gnL 





oOo 
OO 


(b) —— = ames 





— 1¢) | 2 
-- 0 t 2 _ 
(0) PERTURBATION LOADING 


Fic. 5. Panel with shear reinforcement. 


As an example, consider for simplicity a tension panel 
with a cutout interrupting only two stringers and re- 
moving three shear panels as indicated in Fig. 5(a). 
The solution for this panel without reinforcement 
shows that the most highly loaded shear panels are 
(i. 3), t—1,. 1), 1), and (—1, —1). 
that these shear panels are reinforced by the addition of 
plates of thickness ¢’ to the skin of thickness / so that the 
total thickness in these panels is / + ¢’. The perturba 
tion loadings to be placed on the basic structure are 
In Fig. 5(b) are shown the 


Suppose now 


those shown in Fig. 5(a). 
four reinforcing plates or the excess material used for 
reinforcement. These plates are of thickness /’. They 
are loaded with shear flows whose magnitudes are the 
same as the shear perturbation in the corresponding 
panels in the basic structure but with opposite sign. 

The conditions that must be satisfied are: 

(a) The stringer load must be zero at (1, 0), (1, 1), 
(0, 0), and (0, 1) toward the outside of the cutout. 

(b) The shear flow in panels (0, —1) and (0, 1) must 
cancel any distributed perturbation load applied to the 
coaming stringers bounding these panels. (In _ this 
problem there are no distributed perturbation loads 
applied here, and this condition is automatically satis 
fied. ) 

(c) The shear stress in each of the panels (1, 1), 
(1, —1), (—1, 1), and (—1, —1) in the basic structure 
must equal the shear stress in the corresponding rein 


forcing plates shown in Fig. 5(b). 


7 and referring to Tables 3(a) and 3(c) yields, for condition (a), 


— 0.0430¢nL = 0 


Condition (c) is expressed by the equation 


[—0.0921 (Pu LL) a 0.0519( pu 1 —_= 0.0387 (Pu L) + 0.0339 (Pu E) + 0.2017g1 ma 0.080111 +- 


These equations become 


0.0155g¢1, — 0.00949;,] (1/t) = gu(1/t’) 


0.1308hun + 0.6294¢,L = P 


~0.0451 py + [0.12 


79 + (t tl’) |dul = 
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For any given ratio ¢/¢’ this set of equations can be solved for the unknown perturbations );; and g,;.__ The stresses 
due to these perturbations are superposed on the basic stress distribution to give the stress distribution about the 
cutout. 

The analysis is also easily extended to handle the problem of reinforced coaming stringers. It is assumed that 
the stringer load is distributed immediately into the added portion of the coaming stringer such that the stress js 
always given by the force divided by the cross-sectional area. Consider again the tension panel with a cutout 
interrupting two stringers and removing three shear panels. Assume the coaming stringers have reinforcements 
of constant area extending one bay on either side of the cutout. Let the cross-sectional area of the added portion 
be A’ so that the total cross-sectional area of the reinforced portions of the coaming stringers is A + A’. 

In Fig. 6(a) are shown the perturbation loads that are to be placed on the structure. The symbol g;; represents 
the magnitude of the distributed perturbation load applied to that portion of stringer 7 lying between rib 7 and 
ribi + 1. In Fig. 6(b) are shown the loads on the added portions of the coaming stringers which have a cross- 
sectional area A’. The following conditions must hold: 

(a) The stringer load at (1, 1), (1, 0), (0, 0), and (0, 1) toward the outside of the cutout must vanish. 

(b) The shear flow in panels (0, —1) and (0, 1) must cancel any distributed perturbation loading applied to the 
coaming stringers bounding these panels. (There are no perturbation loadings applied here, and this condition 
is automatically satisfied. ) 

(c) The stress in the basic coaming stringer must equal the stress in the excess reinforcing portion at the points 
(1, 2), (0, 2), (1, —1), and (0, —1) and at the following points on the side toward the reinforcement: (2, 2), (—1, 
2), (—1, —1), and (2, —1). 

When condition (c) is satisfied, the excess portions of the coaming stringers (of area A’) can be added to the 
basic structure without violating continuity and can thus replace the perturbation loads on their corresponding 
portions in the basic structure. 

To satisfy condition (a), the following equation must hold [assuming B = 7, see Tables 3(a) and 3(b) |: 


P — 0.5000Pn + 0.3158pn + 0.0534h, — 0.0534 fp — 0.0151 fp» — 0.0713p2 — 0.0266 px 
0.0308q),L —_ 0.06429;,L = 0.0078 Gi a 0.02129 ),L = |) 


Condition (b) is automatically satisfied by symmetry. Condition (c) is satisfied by the following equations: 
At (2, 2) on the side toward the reinforcement 


(P — 0.5000 px» — 0.1514 p2 — 0.0178p2 — 0.0534), —0.0151 py + 0.0713An + 0.0266), + 


0.3991gnuL — 0.1792¢nL + 0.0034¢n.L — 0.01559¢nL) (1A) = p(1/A’) 


and at (1, 2) 


(P — 0.3158px2 — 0.2121 py» — 0.0129p2 — 0.0068px2 + 0.0534pn + 0.0151 pu — 0.3991¢.L — 
0.259391 — 0.0034guL — 0.00999¢,,L)1, A (po + Gul) (1/A’ 


These three equations reduce to 
0.1664 fx» + 0.1308h, + 0.124092 = P 
(0.6692 + (4/A’)|fpo — 0.0294);, — 0.20799, = P 
(0.5476 + (4/A’)]pox — 0.0685py + [0.6718 + (A/A’) ]Gnl = P 


For given values of AA’ these equations can be solved for the unknown perturbations fx», /i,, and gi. The 
stresses due to these perturbations are then superposed on the basic stress distribution to give the stress distri 


bution about the cutout. 


Some calculated results for the reinforced cutouts following is noted: In the case of shear reinforcement, 
just considered are compared in Table B with results although the shear flow in shear panel (1, 1) is greater 
for the same cutout without reinforcement. The rein- in the reinforced structure the shear stress in shear 
forced shear panels were taken to have sheet twice as panel (1, 1) is less, and the shear reinforcement tends to 
thick as the uniform sheet; the reinforced coaming increase the stringer load at station (1, 2). In the case 
stringers were assumed to have twice the cross-sectional of coaming stringer reinforcement, the stringer load at 
area of the uniform stringers. station (1, 2) is greater in the reinforced structure, but 

the stringer stress at (1, 2) is less, and the coaming 

On comparing the results for the reinforced struc- stringer reinforcement tends to increase slightly the 
tures with those for the unreinforced structure, the shear flow in shear panel (1, 1). 
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TABLE B 


ae Stringer Loads Shear Flows- 
Intersection of Rib Stringer Shear Shear 
and Stringer Load Panel Flow 
Unreinforced Cutout 
(1, 2) 1.524P (1, 0) 0 P/L 
(2) 1.168P ch, 2 —0.345P/L 
(1, 4) 1 .OS2P (i —0.046P/L 
(1, 5) 1.049P (1, 3) —0.012P/L 
Cutout with Shear Reinforcement (t’ = t) 
rt. 2) 1.581P (1, 0) 0 Pt, 
(1, 3) 1. 156P § eb —0.514P/L 
(1, 4) 1.072P (1, 3) —0.051P/L 
(1, 5) 1.042P (1, 3) —0.011P/L 
Cutout with Reinforced Coaming Stringers (A' = A) 
(4.2) 1.842P (1, 0) 0 P/iL 
(1, 3) 1.024P ti, 3) —(0).348P/L 
(1, 4) 1. 028P (1, 2) —0.019P/L 
(1, 5) 1.020P (1, 3) 0.000P/L 


Wide Cutouts.—Ilf the cutout extends close to the 
edge of the panel, application of the perturbation solu- 
tions computed for the infinitely wide panel will not 
give accurate results. A simple device is proposed to 
refine the analysis for this problem. The procedure is 
to reflect the panel about its edges and solve the prob- 
lem of an infinitely wide panel with a series of cutouts 
as shown in Fig. 7. This does not increase the number 
of equations in the systems (1) or (2) and (3) but 
merely increases the number of points k, / over which 
the summations must be extended. The shear flows 
in the shear panels just off the edge of the actual panel 
all vanish because of symmetry. 

The number of image cutouts required is dictated 
by the desired accuracy. For most problems it will 
probably be sufficient to consider one image cutout on 
either side of the actual cutout. 


RELATION OF MODEL TO ACTUAL STRUCTURE 


In an actual structure the sheet does carry some di- 
rect stress, contrary to assumption (d) under the sec- 
tion “Basic Assumptions.’’ To account for this, an ef- 
fective area of skin is added to the cross-sectional area 
of each stringer; this effective area usually consists of 
the cross-sectional area of sheet halfway to the adjacent 
stringer on each side of the given stringer, so long as the 
sheet has not buckled. Thus, the stringer area, A, 
to be used in the calculations is taken as the actual 
stringer area plus Df. 

Ordinarily, within the bay or bays where a cutout 
exists, the sheet is removed all the way to the coaming 
stringer. For the portion of the coaming stringer 
within these bays, the effective area of sheet is that 
area extending halfway to the adjacent stringer on one 
side only. Thus there is a discontinuity in the effective 
area of the coaming stringer. This can be taken into 
account in the following manner: Assume that the 
coaming stringer has a reinforcement of negative area 
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Fic. 7. Typical structural model for solution of wide cutout. 


bt/2, and apply the procedure described for analyzing 
cutouts with reinforced coaming stringers. However, 
a simpler but more approximate procedure would be to 
compute the stringer /oads and shear flows neglecting 
this discontinuity in coaming stringer effective area and 
then take it into account when computing stringer 
stresses. 

In the structural model the ribs or transverse stiff- 
eners can be placed where the ribs in the actual struc- 
ture are located. However, if the shear panels are 
very long relative to the spacing between stringers, the 
transverse stiffness of the skin itself acts as continuously 
distributed ribs of finite stiffness and causes the shear 
flow to vary within the panel. When such is the case, 
it is possible that a better approximation to the stress 
distribution will result if fictitious ribs are inserted be 
tween the actual ribs. The ribs in the structural model 
are rigid, thus further investigation is required to de- 
termine just how to locate discrete rigid ribs so as best 
to approximate continuously distributed ribs of finite 
stiffness. Tentatively, it is suggested that it is con- 
venient to have fictitious ribs spaced so that the cutout 
length is an integral multiple of rib spacing. Also, it is 
felt that it is good to have the fictitious rib spacing 
about equal to the stringer spacing so that the shear 
panels are roughly square. 
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Appendix 


THE UNIT PERTURBATION SOLUTIONS 
Derivation of Difference Equation 


The principle of minimum complementary energy will 
be used to derive the homogeneous partial difference 
equations that govern the stringer loads and shear 
flows in a sheet-stringer panel with rigid transverse 
members. Then the loadings of Figs. 2(a) and 2(c) will 
be applied to the structure in the form of boundary con- 
ditions. 

Fig. 8 illustrates the notation for the stringers and 
shear panels. Fig. 9(a) shows the forces and shear 
flows on the portion of stringer j, which is between rib 
zandribi+ 1. The equilibrium of this portion of the 
stringer gives 


(Al) 


ores ities Pi + (Gis = Fe mB = & 


The principle of minimum complementary energy 
states that, of all stress distributions satisfying the 
boundary conditions on stress, the stress distribution 
that exists in an elastic structure is the one that mini- 


energy m*, where 7* = 


mizes the complementary 
(internal energy) — 
acting through the prescribed surface displacements). 
No displacements are prescribed anywhere on the 
structure in this problem so the second term in the com- 


(work done by surface stresses 
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plementary energy vanishes. The transverse stiffeners 
or ribs are assumed rigid. The shear flow is constant 
within any shear panel, and the stringer load varies 
linearly between ribs. For a sheet-stringer panel of 
this type, the internal energy or stress energy is given by 


S.E. = > DY [(L/6AE) (Pas, 2 + PuPiss, » + P,?) 4+ 
i j 
(bL/2Gt)q;;?]_ (A2 


where the integration over the length of the stringer 
between ribs has been carried out. 

The stringer equilibrium of Eq. (Al) must be main- 
tained. Thus, it is required to minimize S.E. with 
respect to P;; and g;; with the constraint relationships 
given by Eq. (Al). This can be accomplished by use 
of the Lagrangian multiplier technique. 
Then the function that is 


Let u;; repre 


sent a Lagrangian multiplier. 
to be minimized with respect to P;; and qj; is 


@=S.E. — > Dd ws (Pus.; — Py + 
i j 
(qi — Qi, j)L] (A8 

Differentiating with respect to P;,; and then with respect 
to g,;and equating the derivatives to zero gives 
08/0P;; = (L/6AE) (2Pi, + Pigs, 5 + Pian, + 

2P ij) bia (hy LJ ~ be) = 0) 
[wil Rte 


These equations, along with Eq. (Al), become 


P41, j + 4P,; - P, Ses 
(GAE/L) (win, ; — my) (A4) 


0®/0qgi; = (6L/2Gt)-2q:; — Mi, 741 L] = 0 


giy = (Gt/b) (wis — wir 541) (A5 
Piss, — Py + (Qiu — Gi, ja)L = 0 (Al 


The p's can be eliminated between Eqs. (A4) and (A5 
to obtain 


(L/GAE) (P41. Poe Vis 1P;, Pee ie P; sg 
Pia, ; — 4Pi; — Pi-s, 3) = (6/Gt) (qi; — Gi-1, ;) (A6 


Now the q’s can be eliminated between Eq. (A6) and 
Eq. (A1). 


6BAiP i; + Ay (Pins, 5 + 4Py + Pin, )) = 0 (A7 


This results in 


where 
B = (E/G) (A/bt) (6?/L?) 


Or the P’s can be eliminated between Eq. (A6) and Eq 
(A1), resulting in 


6BA,;; qij + Ais (iss, j + tq ;; + qd i-1, )=Q (A8 


Eqs. (A7) or (A8) are the homogeneous partial dif- 
ference equations that govern the stringer loads and the 
shear flows in a sheet-stringer panel with rigid ribs. 
The solution to only one of these equations is required, 
for, if the stringer loads are known, the shear flows can 


be obtained by statics and vice versa. Solutions to 
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Eg. (A7) will be obtained in this report. Eq. (AS) 
will not be used. 


Solution for Concentrated Perturbation Load 


One method of obtaining the solution to Eq. (A7) for 
the loading of Fig. 2(a) will be presented. Define the 


function 


h(x) = >> Pi cos jx (A9) 
Multiply Eq. (A9) by cos /x and integrate from 0 to z. 
This gives 


. 


| hx) cos lx dx = ) | P; cos jx cos Ix dx 
0 ) ( 


) 
dv 


(A10) 


The orthogonality relation is 


‘ . T 
cos 7x cos Ix = é,; 
0 s ) 4 


forjand/ ~ 0. Whenj = / = 0, 


| cos jx cos lx dx = [ dx=f 
0 0 


By symmetry, it is seen from Fig. 2(a) that 

Pp. =P 
Therefore, the sum on the right-hand side of Eq. (A10) 
issimply rP;,. A contribution of (4/2)P;, is obtained 


at —/ and, similarly, at +/. At / = O, the sum is the 


full rP;, immediately. Thus, Eq. (A10) becomes 


/ h(x) cos jx dx = xP i; 
7/70 


The function /,(x) is determined in the following man 


(A11) 


ner: Multiply the difference equation, Eq. (A7), by 
cos jx and sum over / from — ~ to ~. 


> 6 BA;:P i; cos jx + 
> Ay (Prat. + 4Py + Piss) cosjx = 0 (Al2) 


Let 


Using the following equations, 


> Pir 41 COS JX = > pi cos (j — 1)x = 


> pi (cos jx cos x + sin jx sin x) 
7 


y oj 008s jx = YD, 
j=- 


piy cos (7 + 1)x = 


j 


> pij(cos jx cos x — sin jx sin x) 


i] 


Eq. (A12) can be written 
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6BA;h (x) + 2(cosx — 1) X 
[hias(x) + 4ho(x) + hya(x)] = 0 


or 


Riasl(x) + Ay(x)-2y(x) + Aia(x) = 0 (A138) 


where 


—3B + 2(cosx — 1) 


(A14) 
3B + (cos x — 1) 


yx) = 
The general solution to Eq. (A13) is 
hy =c(-y + V7? - 1) + a(-y — V7? —- 1) 


Inspection of Eq. (A14) shows that |y| > 1. Further, 
the condition 
lim h(x) = 0 


i—_ 


must hold. So the required solution of Eq. (A13) is, 


for 7 > i, 


h,=al(-y + Vy? — 1)’ 


lA 


for 7 
h,=c(- y- Vy? — 1)’ 


For convenience this solution can be written in the 


equivalent form: for y > 1, 


h, = (—1)' Ce™™ 

for, y < - (A15) 
h,; = Ce~™ 

where 
A = cosh! |4 (A16) 


where the positive value of the inverse hyperbolic cosine 
is taken. 

For all of the values of B in Tables 1—5, the values of 
y(x) are less than —1; thus, the second of Eqs. (A15) is 
used exclusively. 

The coefficient C(x) can be determined from the con- 
ditions at rib 0, the rib at which the load is applied. 
Along rib 0, stringer loads are zero everywhere except 
for 7 = 0. Atj = O the stringer load is p/2. Along 
ribz = 0; then, Eq. (A9) gives 


ho(x) = p/2 


From Eq. (A15), 


ho(x) = p/2 = C(x)e™™° 
So, 
C(x) = p/2, a constant 
Then, 
h(x) = (p/2)e—™* 


and P,, is given by Eq. (A11) as 
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—A(x) 


: p POs aca ; a = cosh A(x) — sinh A(x) = 
Py = e cos jx dx (A17) 
0 


2a —y(x) — V(x) =] 


Using Eq. (A17) can be written in the form 


ae he | . . - . x\y' 
Py = | E + 2 sin? (x/2) 3B + 4 sin? a7 2 sin (x/2) 44 + 3 sin? 2 cos jx dx (AIS 


Solution for Distributed Perturbation Load and fori > 1, 


In order to obtain the solution for a loading of the GBA. Ps + BAP. + 4Pn + Pia) =0 (Ce 
type shown in Fig. 2(c), the effect of this load on Eq. 
(A7) must be examined. From Fig. 9(b), equilibrium 
of the loaded portion of stringer 0 gives 


Define 


Py — Po + (qo — Qo, —)L = qh (A19) k(x) = ps Pj; cos jx (A24 
i] x 
' Eq. (Al) holds in the balance of the panel. Writing we = : 
Eq (A6) fori = 1, 7 = O and diate the q's “th similar to the h,(x) of Eq. (A9). As before, the equation 
C tween the resulting equation and Eq. (A19) yields *™ 
sf, = k(x) cos jx dx (A225) 
6BA,:Pw + A;(P» + 4Pi + Pw) = —GBqL (A20) 0 
Similarly, writing Eq. (A6) for 7 = 0, 7 = 0 and com- is obtained. Using Eq. (A253) as before, the difference 
bining with Eq. (A19) yields equation 
6BA Poo + A;,(Pw + 4 PW 2 P 1.0) = 6BqL (A21) Ri43(x) r ki(x)-2y(x) + k (x) = Q (A26 


It is permissible to consider 7 > 0 because the solu- results. The solution to Eq. (A26) can be written as 
tion for z < 0 is exactly the same as that for 7 > 0 ex- the second of Eqs. (A15) 
cept the sign is reversed throughout the region 7 < 0. 


Eq. (A7) now becomes, fori = 1, k(x) = C'(x)e~*” ' (A27) 
6BAyPs; + Aj (P2; + 4P1, + Po;) = —6BqLéy (A22) where A(x) is given by Eq. (A16). 
The coefficient C’(x) is found by use of Eq. (A22). Multiply Eq. (A22) by cos jx and sum from — © to ©. 
This yields 
7. 6BA;:P;; cos jx + 2(cos x — 1) > (Po; + 4Pi1; + Po;) cos jx = —6BgL 
which can be written as 
GBA; R(x) + 2(cos x — 1) [Ro(x) + 4hi(x~) + Ro(x)] = —6BqL (A28) 
From Fig. 2(c) it is seen that 
Pi; = = Po 
by symmetry. Therefore, from Eq. (A24), 
ki(x) = —ko(x) 
Using this symmetry relation in Eq. (A28) gives 
6B [ko(x) — 3ki(x)] + 2(cos x — 1) [Re(x) + 3hi(x)] = —6BgL 


Substituting from Eq. (A27), the coefficient C’(x) is obtained 


—6BqL 


C(x) = a 
e 


“ je x (6B + 2(cos x —1))+ 6(cos x — 1) =A ISB! 
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Substituting back into Eq. (A25), 


—6BqLe 


cos jx dx 


1 Mo 
Pi; - 
_/ ties (6B + 2(cos x — 1)| + 6 (cosx — 1) — 18B 


which can be transformed to 


3BqL f? 
es or y| 2B ee ee 
- 0 35 + sin* (x/2) 


It is possible to carry out the integration of Eqs. 
(AIS) and (A29) in terms of elementary functions for 
the special case of i = 1 and for all 7. However, the re- 
sulting expressions are not practical for computations 
when 7 is greater than 5 or 6 because they involve small 
differences of large numbers, and, consequently, there 
is a drastic loss in the number of significant figures. 
Integration of these expressions [Eqs. (A18) and 
(A29) | for other values of 2 is difficult. 

These obstacles can be avoided, however. From Eq. 
(A17), for instance, note that the P;,’s for a concentrated 
perturbation can be thought of, except for a constant 
factor, as coefficients of the Fourier cosine series in the 
interval 0 to 27 of the function that takes the values 


” ail (A30) 


between 0 and m and is constructed symmetric about 
7. Similarly, for the distributed perturbation loading, 
the P's are essentially coefficients of the Fourier cosine 
series in the interval 0 to 27 of the function that takes 


the values 


3B + sin? (x/2) + sin (x/2) V 9B + 3 sin? (x/2) 
(A351) 
between 0 and 7 and is symmetric about z. 


The stringer loads in parts (a) and (b) of Tables 1-5 
were computed by means of harmonic analysis of the 


ww of ° 
‘~" cos jx dx 


; ;(A29) 
+ sin (x 2) V9B + 38 sin? (x/2) 


functions (A30) and (A31) based on a collocation pro 
cedure utilizing 24 points in the interval 0 to m (or 48 
points in the interval 0 to 27). This procedure was 
adapted from reference 6. Shear flows were obtained 
from statics with the use of Eqs. (Al) and (A19). Parts 
(c) of Tables 1-5 were obtained by the appropriate 
superposition of the solutions for distributed pertur 


bation loads. 
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Oscillating Airfoils at High Mach Number 


M. J. LIGHTHILL* 
Unwersity of Manchester, England 


SUMMARY 


A simple formula is given for the pressure distribution on an 
oscillating airfoil in two-dimensional flow at high Mach Number 
The formula is expected to be reasonably accurate if the pressure 
on the surface remains within the range 0.2 to 3.5 times the main 
stream pressure. To illustrate the application of the formula, 
some results for symmetrical airfoils performing pitching oscilla- 
tions are obtained and compared with results obtained from exist 
ing theories in the case of high Mach Number. 


NOTATION 


x = coordinate parallel to undisturbed stream 

y = coordinate perpendicular to undisturbed stream, 
measured away from airfoil on each side of it 

y = }(x,t) = equation of airfoil surface (upper or 
lower) at time ¢ 

U = velocity of undisturbed stream 

pi, pi, 41 = pressure, density, velocity of sound, respectively, 
in undisturbed stream 

M = U/a, = Mach Number of undisturbed stream 

w = “effective piston velocity” = (OV /dt) + (UoOY + 
Ox) 

€ = maximum displacement in airfoil oscillation 

w/2r = frequency of airfoil oscillation 

c = airfoil chord 

a = maximum inclination of airfoil surface to the 
stream 

p = pressure at point on airfoil surface 

¥ = ratio of specific heats (assumed constant ) 

y = Y,(x) + Y,(x, t) = equation of ( =) surface 
of symmetrical oscillating airfoils 

Ws, W4 = values of w resulting from Vs, 7, respectively 

Ap = pressure on lower surface minus pressure at same 
point on upper surface 

a(t) = variable angle of attack (in pitching oscillations ) 

m = nose-up moment of aerodynamic forces per unit 
span about axis of pitching oscillations 

Ma, Mg = coefficients defined as (1/p,Ui2c?) (Om/0a)a=a=0 
and (1/p,U\c*) (Om/0@)e=%¢=0 

T = thickness-chord ratio 

Xo = distance of axis of pitching oscillations behind 


leading edge 
h = ratio of Xo to ¢ 
Basic formula: In motions with M/ large and M/[6 + (e€/c) X 
(we/U)] < 1, the surfaces pressures are given to good approxima 
tion by the equation 


w Wy +1) fw\ Hy tijijwy 
Eanes eMC) 4 2E0(0) 
fr ay 4 ay 12 ay 


(1) INTRODUCTION 


A THEORY OF OSCILLATING AIRFOILS in supersonic 
flow is needed, but a completely linearized 
theory has proved too inaccurate. Some second-order 
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terms (part of an ascending series for them, in powers of 
the frequency parameter) have been found by Wylly! 
and by Van Dyke,’ but their methods are complicated, 
and unfortunately their answers do not agree. The 
method of Jones and Skan* takes the thickness of the 
airfoil into account to higher order but not the ampli- 
tude of the oscillations; this method calls for an ex- 
tensive program of computation in any one case. In 
this state of affairs it seems worth pointing out that, 
in the special case of high Mach Number, there is a 
simple theory that should give fairly accurate results 
(better than second-order) even when disturbances are 
rather large. 

In two-dimensional flow past an airfoil at high Mach 
Number, say \/ > 4, the shock waves and expansion 
waves are set at small angles to the undisturbed flow. 
Two consequences of this are worth special notice. 

(a) Gradients in the direction of the undisturbed 
flow are.small compared with gradients perpendicular 
to it. 

(b) Because velocity components parallel to shock 
waves or expansion waves are not changed by them, the 
velocity components perpendicular to the flow are large 
compared with the disturbances to components parallel 
to it. 

Both these facts contribute to the truth of Hayes's 
result‘ that, to a good approximation, any plane slab of 
fluid, initially perpendicular to the undisturbed flow, re 
mains so as tt 1s swept downstream and moves in tts own 
plane under the laws of one-dimensional unsteady motton. 

Goldsworthy’ showed that the relative error in this 
result of Hayes should be of order 1/.f*, both causes 
(a) and (b) above contributing a factor 1/7. Here it 
is assumed that 76 is bounded, say less than 1, where 
6 is the maximum inclination (in radians) of the airfoil 
surface to the stream. This condition, he points out, 
would in most practical applications need to be satis 
fied from considerations of drag. 

Now the result in italics must remain true even if the 
whole flow is unsteady—for example, if the airfoil oscil- 
lates. This is so evident physically from (a) and (b) 
above? that its mathematical deduction from the equa- 
tions of motion, an easy extension of Goldsworthy’s 
work, is here omitted. 

Now as the slab of fluid moves downstream with ap- 
proximately the velocity LU of the undisturbed flow, 


+ Note that it is of little importance to such a slab of fluid, 
while it is being deformed by the wall, whether or not earlier slabs 
have been subjected to the same deformations. 
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OSCILLATING AIRFOILS 


the position of the piece of solid wall bounding it will 
move normal to the stream with a velocity 


w= (OY /ot) + U(OY/Ox) (1) 


where y = Y(x, t) is the equation of the upper or lower 
surface of the airfoil, with y measured away from the 
airfoil in each case. 

Accordingly, the problem of determining the forces 
on the oscillating airfoil at high Mach Number reduces, 
as in the steady case, to the one-dimensional flow prob- 
lem of finding the pressure on a piston moved, with a 
given dependence of velocity w on time, into otherwise 
undisturbed fluid. 

This, by comparison, is a not too complicated prob- 
lem. It becomes particularly simple if a condition like 
Goldsworthy’s is imposed (as indeed it must if the state- 
ment in italics is to be at all accurate)—namely, that 
the magnitude |w| of the piston velocity, Eq. (1), never 
exceeds the speed of sound in the undisturbed fluid. In 
an oscillation with maximum displacement ¢ and fre- 
quency w/27, this condition can be written 


ew + Ub <a, (2) 


where 6, as before, is the maximum inclination of the air- 
foil surface to the stream. This condition could be 


rewritten 
M[6 + (e€/c) (we/U)] < 1 (3) 


where J = UL//a,; is the Mach Number and c is the 
chord, so that we/ Ll’ is the frequency parameter. Evi- 
dently, condition (3) might well be satisfied in practical 
problems. (The theory would still have value as a 
rough approximation even if the left-hand side some- 
what exceeded 1.) As will be seen below, it per- 
mits compression up to three and one-half times and 
rarefaction down to one-fifth of the undisturbed pres- 
sure. 

The advantage of imposing this condition is that the 
pressure on the piston then depends, with fair accu- 
racy, only on its instantaneous velocity. The pressure 
distribution is in fact a simple function of the piston 
velocity w, given by Eq. (1). 

When there is no shock wave, indeed, this is true 
without any restriction on amplitude, and the pressure is 


given by the “‘simple wave’’ condition 


— Ly\a—! 
p = (; + L ) (4) 
pi 2 a 


where 7 is the ratio of the specific heats (assumed con- 
stant). When there is an initial shock wave (because 
the piston moves initially with finite velocity into the 
fluid), Eq. (4) can be improved by use of the ‘‘shock- 
expansion” theory, which takes into account the exact 
pressure change at the shock [which is less than Eq. 
(4)] and assumes a simple wave (i.e., neglects lateral 
entropy gradients) behind the shock. But this theory 
gives pressures depending not only on the instantaneous 
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TABLE | 
w/a, Simple Wave Shock Expansion Quadratic Cubic 
1 3.583 3.473 3.240 3.520 
0.0 1.949 1.916 1.910 1.945 
0 1.000 1.000 1.000 1.000 
—0.5 0.478 0.488 0.510 0.475 
0.220 0.440 0.160 


—] 0.210 


value of w/a, but also on the value that w/a, had when 
the shock wave was produced. 

In Table | is shown the dependence on w/a, (for y = 
1.4) of (a) the simple wave expression, Eq. (4), for 
b/pi; (b) the value of p/p; on ‘“‘shock-expansion”’ 
theory, taking the shock at its maximum permitted 
strength, corresponding to w/a, = 1; (c) the “Buse- 
mann’ quadratic approximation—namely, the first 
three terms of the binomial expansion of Eq. (4); (d) 
the cubic approximation (first four terms of the ex- 


pansion )-—namely, 


w (y + 1) /w\? 
otter eo (“) + 


Pi ay | ay 
v(y¥ + 2 (2) 5 
(.>) 
12 ay 


Since for a given value of w/a, the value of p/p, may 
vary continuously between the extremes of the “‘simple 
wave” and “‘shock-expansion’’ values (as the value of 
w/a, at the shock varies from 0 to 1), any approximate 
value that always lies reasonably close to both is ac 
ceptable. The cubic approximation, for example, 
which is always within 0.06 p, of both values (not 
a great error where such large pressure differences are 
occurring), even if w/a; is as high as 1, is in many ways 
to be preferred even to the exact simple wave value, 
which is in error by 0.11 /; immediately behind the 
Since, also, the cubic is reasonably 
(Note that cubic 


shock wave. 
simple, it will here be adopted. 
terms in the expansion of pressure ratio in a simple 
wave have often been discarded on the grounds that 
their coefficient is not the same as in the expansion of 
the pressure ratio through a shock; however, the co- 
efficient in the latter case 1s y(y +. 1)°/32, which 
differs from that in Eq. (5) by only 10 per cent for + 
1.4, so that the terms are well worth retaining. ) 


Finally, then, the pressure at any point on an oscil 
lating airfoil will be taken as given by the cubic (5) in 
w/a,, where w is defined in terms of the equation y = 
V(x, ¢) of the airfoil surface by means of Eq. (1).* 

It is unnecessary for the author to make extensive 
calculations here on the basis of the postulated pres- 
sure distribution, since it is a simple matter for anyone 
concerned in technical applications to make use of it to 
obtain any particular results he may need on the pres 


* It should be noted that a formula of somewhat similar char 
acter was suggested by Jones® as an empirical approximation for 
general M. The present work indicates that Jones’s formula is 
right for large M, and so gives it some added support 
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sure forces acting on oscillating airfoils under the con- flower surface| 


ate ; : é y= Vis) = Viet) - : (6) 
ditions here treated. Some specimen results will, (upper surface { 
however, be displayed and compared with existing 
theories. Hence, by Eq. (1), 

(2) SYMMETRICAL AIRFOILS kt oy; _OY,; 
w= [ } . (¢) = + l _ 

Supersonic airfoils are usually symmetrical, so that or Ox 
at zero angle of attack the equation of both surfaces w(x) = w(x, t) (7) 


[with y measured away from the airfoil in each case, 

as in Section (1)] is the same, say y = Y,(x). Changes Now in many problems one would wish to know not 
in y due to oscillations, whether pitching or bending, the pressure distribution, Eq. (5), on each surface but 
will have opposite sign on both sides, say Y;(x, ¢) on only the load distribution, obtained by subtracting the 
the lower surface and — );(x, ¢) on the upper surface. pressure distribution on the upper surface from that on 


The full equation of both surfaces is then the lower surface.. This is 
Ap w, + Ww; vy (y + 1) fw, + wi\? y (vy + 1) (ws + wi\? 
=1+y¥ a ‘ fe es 
Pi ay 4 ay 12 ay, 
Wy, — W; y¥(y¥ + 1) (* = my" y¥ (y+ 1) (*: ~ my 
—% pa = : 
ay f ay 12 ay 
Ws y(y + 1) /w,\* ] w, l w,\° ’ 
= fort tn + “ ( )] +-+l9 + (8) 
ay 2 ay ay 18) ay 


As would be expected, there are no terms in the square of the oscillatory part w; of w. The main effect of nonlin- 
earity is to change the coefficient of w; a, by a factor depending on Mach Number and on the local slope of the 
airfoil contour. 

It will be sufficient as an example to consider a pitching oscillation, say one about the axis x = xo, vy = 0, with 
variable angle of attack a(t). Then 


Y,(x, t) = a(t) (x — xX), W, = a(x — x) + Ua (9) 


Substituting also for w, in Eq. (8), it becomes 


Ap - Vly + 1) cece re a(x — Xo) 
= 127+ y¥(7¥ + 1)MY,'(x) + : M?Y,’?(x) + Ma}+ 
d ay 


Pi 
¥ ( 1) | a(x — xo) 8 
nyt k “+ Ma| (10) 
6 a; 


The terms in (10) which are linear in a and & are especially important for determining the aerodynamic stiffness 
and damping, respectively, of the oscillation. The stiffness is the coefficient of a, and the damping is that of a, in 


—m = / (x — x) Ap dx (11) 
0 


of the load distribution about the axis x = x). (Here, c is the chord.) 
Thus, by Eq. (10), the aerodynamic stiffness coefficient is 


l on) | if | y¥+1 
ira = : = = (x — x) |] 2 (4 1) MY,’(x) + M?Y,’2(x) | dx (12) 
pUj"c? ( OES cacimd Me? y| . . ry) “ 2 1 it 


and the aerodynamic damping coefficient is 


l om l “ y+ 1 : 
—m, = —{— = 2 (x — x9)? } 2 (4 1) MY,"(x) + M?Y,'*(x) f dx (13) 
pil =f .. ; aa | la: | 5 ' 2 at 


The latter is essentially positive—that is, there is always damping at high Mach Number. (This is obvious from 
the point of view of this paper, the damping being essentially the radiation damping of a piston in air.) 
To fix the ideas still further, consider a biconvex airfoil of thickness-chord ratio 7, given by the equation 


the nose-down moment, 


V(x) = 2r[x — (x?/c)] (14) 
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(Note that, by the restriction (3), we must have 7 < (2.1/)~' even for small oscillations; somewhat less latitude 


can be allowed if the oscillations are large.) Then writing x)/c = h for the position of the axis of pitching oscilla- 
tions relative to the chord, Eqs. (12) and (13) become ‘ 
7 I a 2.2 2 — 
—m, = s1—-(y t+ LMe + = (4 + ODM E?§ — hI 2 + - (+ DM? Ie (15) 
M\ 3 3 3 | 
1 42 


l 4 l l 
~~ (9 % 1)Mr+— (vy + 1)AMPr? - on ~-—ie + fees — 9 + 1a | + 


“me V3. 3 15 


») 


h? E + — (y + Waal (16) 


The value of 4 for which Eq. (15) vanishes, corresponding to the aerodynamic center, is tabulated against Mr 
for y = 1.4in Table 2. This table also gives the value of —m, in the special case (hk = 0.5) of oscillations about 


an axis at mid-chord, when 


l l | 
—m, = ( 1) M?,7? 7 
. gle+. 1+ ") (4) 


Now the expression given by Wylly' for this last quantity is, in our notation, 


—m, = 


M? —2 
eet sa : ,[ sar + (2y — 16)M* + (7 + 19)A22 - 4] (18) 
WN _ Wu ae kh 


This gives, as .\J — o for fixed Mr, 
—m, ~ (1/M) [(1/6) + (2/3) Mr] (19) 


Now Wylly does not retain terms of high enough order to give a term in \/*r? in — \/m,, but the coefficient of Mr 
ought to be correct. However, it is seen from Eq. (17) that, in this case of h = 0.5, the coefficient of Mr in — Mm, 
for large .\J should be zero. 

The theory of Section (1), if no slip has been made (and it is very simple), must give correctly the asymptotic 
form as .\J —~ © of any term of second order in the disturbances, since there is no difference between the simple 
wave and shock-expansion theories to this order. (In contrast, the third-order terms, which have been retained, 
should merely improve the accuracy, without being mathematically exact.) Therefore, the conclusion seems in- 
escapable that some error has crept into Wylly’s calculations. 

Van Dyke's expression’ for the pressure on the upper surface may be written, in our notation, as 


en a d = 4 y ap) 2 — Yep op 
p ats ae é =") 4. (1 t)y(¥ ee = yM + Y (* mnt 4 
) Ay? — 1)2 


pi (M? — 

6 M M?(2 + (vy — 1)M | 
=e - 7,85 | a : - | E V(x) + - xMY,’(x) l (20) 
a, \((M? — 1)” (M2 — 1)3 2 f 


a,;? 


To obtain this result, he neglects w;*? and all the third- (from 0.062 to 0.072 as the frequency parameter varies 
order terms and takes the limit as the frequency param- from 0 to 0.75), do not differ enormously from Eq. (17), 
eter tends to zero. For large M the last term in Eq. which has the value 0.084 (independent of frequency 
(20) vanishes, and the rest becomes identical with the parameter). Since the Jones and Skan values for —m, 
basic formula of this paper (which holds independently are expected to have good accuracy, this result at a 
of frequency parameter) if, in it, all second-order terms moderate Mach Number augurs well for the present 
except w,”, but no third-order terms, are included. This 

constitutes a valuable check on Van Dyke’s work. TABLE 2 


The vreates ac Number fo yhic ones and , - 0 . ; ; : 
Phe greatest Mach Number for which J . . Distance of Aerodynamic Center from Leading Edge as Fraction 


Skan* have carried out their calculations is JJ = 2, of Chord, and Aerodynamic Damping Coefficient for Pitching 


: . Oscillations About a Mid-Chord Axis, as Functions of the Product 
at whic agreeme y > prese heory would : a “ 
t which poor agreement with the present t ’ of Mach Number and Thickness-Chord Ratio for a Biconvex 


be expected. However, their results do agree with it in Airfoil at High Mach Number 
predicting that the dependence of —m, on the thick- yy, 0 0.1 0.2 0.3 0.4 05 
ness-chord ratio 7 for oscillations of a biconvex airfoil {hk}, 0.500 0.460 0.422 0.388 0.358 0.333 
about a mid-chord axis becomes negligible (for moder- Eee Eee 0.167 0.169 0.176 0.188 0.205 0.227 
iat VU M WU VW M M 


ate 7) at lJ = 2. Indeed, their actual values of —m, _ 
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theory at the rather high Mach Numbers for which it REFERENCES 
was designed to work. 


in a and & in Eq. (10), not so far discussed, the one that California, 1951. 
will predominate when the frequency parameter is less ? Van Dyke, M 
than | is 


., Private communication. 


3 Jones, W. P., and Skan, S. W., Aerodynamic Forces on Bi. 
( 6)-y(y + 1)M%a3 (21) convex Aerofoils Oscillating in a Supersonic Airstream, ARC 
Reports 13162 and 14217 (limited circulation), to be published as 
This uniform load distribution along the airfoil consti- A.R.C. R. & M. 2749. 
tutes a nonlinearity in the stiffness for a pitching os- * Hayes, W. D., Q. App. Math., Vol. 5, pp. 105-106, 1947 
cillation, which will help to stabilize it when the axis is § Goldsworthy, F. A., Q.J.M.A.M., Vol. 5, pp. 54-63, 1952. 


ahead of the mid-chord position (but may render pos- 6 Jones, W. P., The Influence of Thickness-Chord Ratio on Super. 


sonic Derivatives for Oscillating Aerofoils, A.R.C. Report 1087] 
it for axes to the rear of it. (limited circulation ), to be published as A.R.C. R. & M. 2679. 


sible some subharmonic resonance) and to destabilize 





i) 
\ — 7 a ‘. . ° ° 
Che Upwash Correction for an Oscillating 
; : ter on 
, Wing ma W ind | unnel 
} (Concluded from page 386) 
7 everywhere except in a small region about g = k; 
J0 here f(g) may be replaced by f(k) and brought outside 
If w appears in the numerator and (g + &) in the the integral, provided f(g) is continuous. The result- 
denominator, then, as u — 0, the integral vanishes. ing integral may be evaluated, and then, as u > 0 
If « appears in the numerator and g — k in the the result is rf(k). This is a special case of Fourier’s 


denominator, then, as uh — O, the integral vanishes integral theorem. 





: , . ‘Wylly, A., A Second-Order Solution for an Oscillating, Typ. 
One may note, in conclusion, that of the cubic terms Dimensional, Supersonic Airfoil, The RAND Corporation Report 
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The Contribution of Higher Mode Resonance 
to Helicopter Rotor-Blade Bending’ 


HAROLD HIRSCH? 
Hughes Aircraft Company 


ABSTRACT 


Results obtained from flight tests of a family of main rotor 
blades, similar in all respects except flapwise bending stiffness, 
are reviewed with the objective of pointing out the role that res- 
onant bending, particularly in the second and third bending 
modes, has been observed to play in determining the beam 
bending moments experienced by the test blades. It has been 
found that in certain critical conditions the contribution of 
higher mode bending to the total bending moment may be as 
great as the sum of the steady, first and second harmonic terms 
usually considered in current blade analysis methods. The source 
of the forces exciting the resonance has not been determined; 
however, it appears to be of aerodynamic origin. 

The results call attention to the influence resonant bending 
may have in determining, from the fatigue point of view, the 
structural lite of rotor blades. Since the blade resonant fre 
quency in the second mode will, in general, lie between four and 
five times rotor speed and in the third mode between six and 
nine times rotor speed, the number of stress cycles experienced 
by the rotor blades will include corresponding multiples of the 


rotor rotation 


INTRODUCTION 


ee THE GREAT AMOUNT OF EFFORT already de- 
voted to their study, there are many important 
helicopter problems that are still lacking final solution. 
Included among these are the interdependent problems 
of rotor dynamics, rotor-blade structural analysis, and 
helicopter stability and control. By applying mathe- 
matical methods and physical principles developed in 
the airplane and other fields, considerable progress has 
been made toward their solution. The principal short- 
coming of methods for treating these problemis is, in the 
present state of the art, a lack of adequate procedures 
important aerodynamic loads 
The practical problems that 


for determining all 
applied to the rotor. 
arise when attempting to predict realistically the aero- 
dynamic loads result principally from two sources: 
the complexity of induced flow effects and the inter- 
dependent manner in which rotor-blade motions and 
structural deflections enter into the determination of 
the aerodynamic loads applied to the blade. 
Recognizing the important role played by blade 
bending, several investigators have undertaken to 
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Rotary Wing Branch, Propeller Laboratory, Wright Air Develop- 
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t Work done while at Cornell Aeronautical Laboratory, Inc. 
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introduce the effects of blade bending deflection into 
the analytical treatment of helicopter problems. 
Johnson and Mayne! have outlined two methods for 
introducing the inertial effects resulting from blade 
bending deflection into blade structural analyses. 
One method solves the problem by means of a lengthy 
“tabular-dynamic’’ method; the other obtains sub- 
stantially the same results with much less work by 
multiplying the ‘‘static’’ blade moments by the ampli- 
fication factor for a simple undamped single degree of 
freedom system. Winson? has proposed a method 
applicable to rigid rotor blades which takes into ac- 
count the effect of first- and second-mode blade bend- 
ing deflections in determining the blade aerodynamic 
loads. Flax and Goland* have described both tabular 
and amplification factor methods for introducing the 
additional inertia terms and the aerodynamic damping 
effects due to blade bending deflections into blade 
structural analysis procedures. Finally, Arnold and 
Goland* have developed a method that introduces the 
effects of blade flexibility into helicopter stability and 


control analyses. 


Attention has been called to the foregoing references 
in order to show that awareness of the potential role 
that might be played by rotor-blade bending, particu- 
larly in the first bending mode, has existed for some 
time. However, the conclusions derived from these 
methods have been only as realistic as the assumptions 
on which they are based. In practice, this reduces to 
the fact that if knowledge of the important bending 
modes and their associated aerodynamic effects is not 
and that the role of these factors is, there- 
theory 


available 
fore, not included in the formulation of theory 
will never, of itself, predict their existence. It is 
important, accordingly, that experimental information 
be made available which will show which effects are of 
practical importance and should be included in theo- 
retical formulation. Further, if important phenomena 
are currently being overlooked because of the fact that 
they are not indicated by existing analytical methods, 
it is even more important that experimental programs 
be undertaken so that their existence will be discovered. 


The research program from which the subject matter 
of this paper has been extracted was undertaken with 
the general objective of providing comprehensive 
experimental information that would satisfy the dis- 
cussed need for factual data concerning rotors. Spe- 
cific objectives of the program were to determine the 
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Fic. 1. The three sets of fiberglass skin and balsa-wood core 
blades flown in the program. Stgain gages are wired through 
plastic tubes built into the blades. 





BLADE NODE LOCATION 
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Fic. 2. Diagram showing node positions for each blade and 


giving blade spanwise stations for bending data channels. Sta 
tions are measured from the center of rotation. 


bending moments and motions experienced in flight by 
articulated rotor blades and thereby to explore possible 
causes of a number of unexplained rotor-blade fatigue 
failures. It is believed that the program was highly 
successful in realizing its objectives, particularly since 
the fundamentally important rotor-blade resonance 
phenomenon described in this paper was disclosed in 
the course of the project. It is believed that the re- 
sults obtained, in addition to clearly revealing the physi- 
cal nature of the resonance problem, call attention to 
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the need for re-evaluation, or extension, of current cop. 
cepts of rotor dynamics and aerodynamics. Although 
the source of the forces exciting the measured reso. 
nances is not yet known, it is hoped that additional 
investigations will soon provide a satisfactory answer 
to this problem. 


BACKGROUND 


Program Initiation 


The so-called ““H-5 Variable Stiffness Blade Pro- 
gram,” from which the results presented in this paper 
are extracted, was conceived and planned jointly by 
W.A.D.C. Rotary Wing Branch personnel and engi- 
neers of the Cornell Aeronautical Laboratory. Ap.- 
proaching the problem with the limited theoretical and 
experimental data then available, it is understandable 
that the existence of higher mode resonant bending was 
not foreseen at the time the program was undertaken, 
In fact, at the time there were no proved arguments to 
refute the popular opinion that the role of structural 
blade bending stiffness in determining blade motions 
and stresses was of secondary importance. This 
opinion, while now proved false as a generalization, can 
be justified with some logic for the case of first-mode 
bending by the fact that the centrifugal stiffening 
effect realized with normal blades will be four or five 
times as great as that due to their structural stiffness. 
In view of this fact and anticipating only quasi-static 
bending in the first bending mode, with some ampli- 
fication, it was felt that the differences between data 
measured with blades possessing the limited amount of 
stiffness variation obtainable with conventional types of 
blade construction would be too small to be determined 
accurately. Design requirements were, therefore, 
established for blades that would provide a four-to-one 
stiffness range and still possess approximately the same 
mass and aerodynamic characteristics. Further, the 
blades were to possess one-half, one, and two times the 
stiffness of standard H-5 blades and be flyable on an 
H-5 helicopter. With a four-to-one stiffness range, it 
was anticipated that stiffness variation effects would be 
sufficiently large so that they could be measured. 


Fiber-Glass Blade Design 


Extensive preliminary studies’ of various blade de- 
signs produced only one that met the four-to-one stiff- 
ness range requirement. The proposed design called 
for the development of three separate sets of blades, one 
of each stiffness, employing a fiber-glass shell, die 
molded integrally to a balsa-wood core. By employing 
the high strength-to-weight ratio of fiber glass, in con- 
junction with the ability to vary the elastic modulus of 
the fiber-glass skin by selectively orienting the fiber 
direction, the required stiffness range was provided. 
Whirl test blades were first constructed to prove the 
fiber-glass design, and, after they had _ successfully 
passed Propeller Laboratory tests, the three sets of 
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flight-test blades were fabricated. Detailed informa- 
tion concerning the design and manufacture of the fiber- 
glass blades is given in references 6 and 7. A further 
description of the blades, together with an analysis of 
their mass, structural, and dynamic characteristics, is 
also available. The finished appearance of the 
blades is shown in Fig. 1. 


Flight-Test Program 


In planning the flight program, an attempt was 
made to select a range of test variables sufficiently 
large so that the conditions investigated would be 
representative of both normal and abnormal operation. 
The steady-state program that was flown is shown by 
Table 1. Tests with the low-stiffness (No. 1) blades 
were not carried to high advance ratios or into climb 
and autorotation conditions. These tests were elimi- 
nated for reasons of safety, since resonant bending 
characteristics of the low-stiffness blades were particu- 
larly pronounced. 

One blade of each stiffness was instrumented to 
measure bending at nine spanwise stations. The sta- 
tions selected are shown in Fig. 2, which also shows 
experimentally determined stationary node locations 
for first-, second-, and third-mode bending. In addi- 
tion to the beam bending traces, the test records con- 
tain traces for the pitch, flap, and lag motions of the 
instrumented blade. The load developed in the lag 
damper of the instrumented blade and the load due to 
blade pitching moment developed in the vertical con- 
trol rod to each of the blades was also recorded. Only 
the blade bending data are considered in this paper. 


Instrumentation 

Since the variable stiffness program was largely con- 
cerned with interpretation of the differences between 
measured quantities rather than the absolute magni- 
tude of the quantities, an important and effort-con- 
suming part of the program consisted of the design, 
development, and accuracy investigations required to 
provide adequate instrumentation. Realizing that the 
program results would only be as reliable as the instru- 
mentation, studies were undertaken for the purpose of 
evaluating alternative instrumentation possibilities. 
In these studies, emphasis was placed on devising a 
system that could simultaneously record as many as 
28 channels of dynamic data measured on the rotating 
components of the rotor system. Results of the 
studies indicated that a multiplexing system — with the 
subcarrier oscillators, bridges, and amplifiers mounted 
on the rotating head—would satisfy the program re- 
quirements. In the completed system, the measured 
data, after being amplified and multiplexed, were trans- 
mitted to a ground-located receiving station where 
they were fed into analyzer circuits that resolved the 
combined signal into the individual data channels. 
Each channel was passed through a low-pass electrical 
filter to reduce hash and then was recorded on an 
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WAN ROTOR HEAD INSTRUMENTATION 


The 28-channel rotating telemetering dome assembly 
and blade root data pickups. 


Fic. 3. 


oscillograph. A view of the rotor head, telemetering 
dome assembly, and several data pickups is shown in 
Fig. 3. Detailed discussion of the design, operation, and 
accuracy of the instrumentation can be found in refer- 
ences § and 9. 

The instrumentation 
employed in the program, contributing the following 


developed was_ successfully 


advantages: 
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(1) Substantial reduction in weight and size of air- 
borne equipment. 

(2) Elimination of unamplified signal transmission 
through slip rings. 

(3) Reduction in number of slip rings required per 
channel. 

(4) Simultaneous recording of a large number of 
channels. 

(5) Continuous telemetering of flight information to 
the ground station, where engineers could observe the 
data that were being recorded. Malfunctioning equip- 
ent, off-scale traces, loss of signal, and excessive 
strains were often detected by ground observers while 
the flight was in progress, and corrective action was 
taken before data were recorded. 

Disadvantages of the system consisted primarily of 
the need for excessive maintenance due to lack of 
ruggedness and Both of 
these shortcomings are attributed to the fact that the 


temperature instability. 
program did not allow for complete development of 
the system, rather than to any inherent system charac- 
teristic. In the program, correction for temperature 
drift was introduced in the data reduction procedure, 
using correction factors obtained from “‘calibrate inject” 
signals recorded before, during, and after each flight. 
Temperature instability was, however, a serious prob- 
lem with the instrumentation. Experience based on 
this and other prograins indicates that too much atten- 
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Reproduction of blade shakedown oscillogram for the medium-stiffness blade from which higher (third-) mode resonant 


bending was first identified 
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tion cannot be given to anticipating and eliminating 
temperature effects when designing flight-test instru- 


mentation. 


Discovery of Resonance Phenomenon 


Prior to installation of the test instrumentation, each 
set of blades was subjected to ground and flight shake- 
down tests. Lacking instrumentation, these prelimi- 
nary tests were evaluated on the basis of pilot’s opinion. 
It was reported that the medium-stiffness (No. 2) 
blades appeared to possess normal characteristics. 
The high-stiffness (No. 4) blades were reported to be 
more responsive to the controls but, at times, somewhat 
“rough.” The low-stiffness (No. 1) blades were judged 
to be slow to respond to control inputs and to be gener- 
ally rough, often producing a cyclic stick “chatter’’ 
during the application of control motions and occasion- 
ally producing cyclic stick shake at high frequencies. 

Being concerned with these observations concerning 
flight characteristics of the blades, it was decided that 
additional blade shakedown tests would be conducted 
when installation of the instrumentation was com- 
pleted. Instrumentation was subsequently completed, 
and a ground run-up was made with the medium-stiff- 
ness blades. A section of the blade bending records 
obtained is reproduced in Fig. 4. 

Glancing over the run-up record, it was observed 
that a periodic oscillation, prominent in the figure, was 
found to appear and disappear throughout the records. 
Noting the regularity of the oscillation and its uniform 
appearance, it was presumed that the instrumentation 
was malfunctioning. Since the oscillation appeared to 
be common to all bending channels, it was anticipated 
that the difficulty would be located in some part of the 
instrumentation carrying the multiplexed signal, such 
as the master amplifier or the radio links. A check 
showed, however, that all instrumentation was operat- 
ing properly. 

The oscillograph records were then more critically 
re-examined, revealing that in reality only nine of the 14 
recording channels were oscillating. This might have 
been rationalized on the basis that the sensitivities of 
the remaining channels were too low to show the 
oscillation. It was also discovered, however, that at 
least one of the channels (No. 9) was unquestionably out 
of phase with the others, and it appeared that channel 
No. 10 was similarly out of phase. 

Becoming aware that something unanticipated was 
taking place, an intuitive step was taken. Picking an 
instant, indicated by the vertical line drawn on the 
oscillogram, when the periodic trace displacements were 
maximum, the displacement direction of each trace was 
tabulated. It was noted that the higher order content 
of trace No. 2 was small and that trace No. 8 appeared 
to be relatively free of higher order components. It 
was also observed that the higher order component in 
all cases appeared to be eighth order. 
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Fic. 5. Resonance diagram for flapwise bending of the fiber 


glass blades in the first, second, and third bending modes 


Next, as shown in the figure, a base line was drawn 
and the relative spanwise location of each trace was 
approximated. The direction of trace displacement 
was then indicated by an arrow, letting traces No. 2 and 
No. 8 remain undetermined. Knowing that all chan- 
nels were sensed in the same direction, the blade deflec- 
tion pattern immediately identified the source of the 
eighth-order component; the blade was unquestionably 


bending in the third bending mode. 


Having thus revealed the existence of higher (third) 
mode resonant bending, steps were taken to determine 
whether its occurrence could be correlated with the 
measured frequency. Using available methods’ " 
this was readily done, as shown in Fig. 5, and later dis- 


cussed. 


Natural frequency calculations were subsequently 
undertaken for the low- and high-stiffness blades for 
the cases of first-, second-, and third-mode bending. 
The results obtained predicted that at least one par- 
ticular combination of blade stiffness, mode shape, and 
rotor speed would produce a resonant condition corre- 
sponding to each harmonic order between two and nine. 
Thus, before the flight program was undertaken, the 
resonant conditions and frequencies were predicted in 
detail. 
paper will show how well the predictions were checked 


Results presented in later sections of this 


by the measured data. 








412 TOURNAL OF 


THE AERONAUTICAL SCIENCES 


JUNE, 1953 


TABLE 2 
Summary of Nonrotating Blade Characteristics 


Blade Stiffness Weight, Lbs. 


Low (No. 1) 162.1 
Medium (No. 2) 161.7 
High (No. 4) 160.7 


Ist bending, 2nd bending, 


c.p.m. c.p.m. 
155 440 
200 535 
260 660 


DETERMINATION OF BLADE NATURAL BENDING MODES 
AND FREQUENCIES 


As the first step in obtaining a correlation between 
the resonance conditions measured and those analyti- 
cally predicted, the natural blade beam bending fre- 
quencies were determined. As previously mentioned, 
this was done for the first three bending modes for each 
blade stiffness, using both Horvay’s” approximate 
coefhicient method and the tabular solution outlined by 


Myklestad."! 


Horvay Solution 


In obtaining an approximate solution with Horvay’s 
rotational coefficients, experimentally determined sta- 
tionary blade natural bending frequencies were em- 
ployed rather than those computed from Horvay’s sta- 
tionary coefficients. The experimentally obtained fre- 
quencies were used since it was believed that they 
would be more correct than those computed from the 
coefficients. This point, however, is not too important 
since the natural frequency values in the rotor operating 
speed range are insensitive to moderate variations of the 
stationary frequency values. This is particularly true 
of the first-mode frequencies and less so for the second- 
and third-mode frequencies. 


Experimental Determination of Stationary Natural Fre- 

quencies 

The stationary frequencies used in the Horvay solu- 
tion were determined from a vibration survey of the 
blades. Substituting a welded steel fitting for the 
standard H-5 root sleeve and flapping link fittings, the 
blades were cantilevered from a structural backstop 
with the spanwise feathering axis horizontal and the 
chord plane vertical, leading edge down. The welded 
fitting also served to extend the blade so that the flap- 
ping hinge was located at its normal position at Station 
0. The blade was hinged with a vertical pin at the 
attachment to the backstop so that it was permitted 
freedom in the flapping plane while being rigidly re- 
strained in torsion and against movement in the blade 


lag plane. A variable-speed mechanical vibrator was 


Moment of Inertia 


Spanwise About Station 0, 


C.G. Station, In. Slug-Ft.? 
123.3 740 
123.5 745 
125.0 747 


—Blade Natural Frequency 


3rd_ bending, Ist torsion, 


¢c.p.m. ¢c.p.m. 

895 1000 
1110 1110 
1290 1310 


attached to the blade tip at the quarter chord, and fre- 
quency sweeps were made to determine the frequencies 
of the first three natural bending modes and the node 
positions for each. The results obtained from the 
blade vibration survey are presented in Fig. 3 and 
Table 2. 


Measured weight and balance characteristics of the 
blades are also summarized in Table 2. It will be noted 
that the maximum weight deviation of any instru- 
mented blade from the average value for the three 
blades is less than | per cent and that the maximum 
spanwise center of gravity deviation from the average is 
More important, the varia- 
Since the 


also less than | per cent. 
tion of moment of inertia is even smaller. 
measured node locations for all blades (Fig. 3) are in 
good agreement and since the type of blade construc- 
tion precludes large local variations in blade weight, it 
may be concluded that the spanwise distribution of 
bending stiffness (//J) is similar in all blades, even 
though they differ widely in bending stiffness magni- 
tude. 


Myklestad Solution 


The Myklestad method consists of a tabular proce- 
dure for determining the frequency at which the blade 
will vibrate subject to inertia forces while satisfying 
boundary conditions of zero deflection and zero moment 
at the hinged root and an assumed unit deflection at the 
blade tip. Trial frequencies other than a natural fre- 
quency produce a residual moment at the root hinge. 
By plotting the residual moment at the root hinge ver- 
sus trial frequency, a curve is obtained from which the 
frequency at which the residual moment for each mode 
This frequency is 
The use 


is zero can be accurately selected. 
the natural bending frequency for the mode. 
of the preliminary natural frequency values deter- 
mined from the Horvay solution made it possible to 
choose the trial values of frequency near the natural 
values, thereby reducing the number of trial computa- 
tions required. This procedure is recommended as a 
timesaving step. 





Con 


que 
My 
Fig 
tair 
rest 
mo 
gen 
The 
cur 
is 1 
rot 
906 


err 


pro 
be 

put 
spr 
sha 


Cor 


giv 
Aer 
que 
nat 
equ 


wh 


sul 
giv 
cal 
Sta 
res 
col 


| fre- 


ncies 
node 
the 
and 


the 
oted 
‘tru- 
hree 
lum 
re 1S 
ria- 
the 
2 in 
ruc: 
E at 

of 
ven 
ni- 


Ce- 
ide 
ing 
nt 
he 
re- 


er- 








RESONANT BENDING AND H 


Comparison of Blade Frequency Results 


The relationships between blade natural bending fre- 
quency and rotor speed obtained by the Horvay and the 
Myklestad solutions may be compared by reference to 
Fig. 5, which includes curves showing the results ob- 
tained by each method. As indicated on the figure, 
results are given for the first, second, and third bending 
modes for each blade stiffness. It is seen that, in 


general, the agreement between methods is good. 
The maximum error between any two corresponding 
curves is less than 5 per cent for the stationary case and 
is less than 3 per cent in the range corresponding to 
rotor speeds at which flight tests were conducted (175 
206 r.p.m. It is also observed that the direction of 


error appears to be random. 


Since the difference between results is less than the 
probable overall accuracy of either calculation, it will 
be assumed in subsequent discussions that the com- 
puted natural frequencies lie somewhere between the 
spread of the two results——that is, they lie within the 
shaded areas between curves. 


Comparison of Blade Frequency Coefficients 

Blade frequency coefficients comparable to those 
given by Horvay have been calculated for the Cornell 
Aeronautical Laboratory blades, using the natural fre- 
quencies obtained from the Myklestad solutions. The 
natural bending frequency may be approximated by the 
equation!” 


Wn? = Wsn? + b, 2? cycles per sec. 
where 


w, = rotating natural frequency in the mth mode 
(cycles per sec. ) 

Ws, = Stationary natural frequency in the mth mode 
(cycles per sec.) 

b, = a constant for the nth mode 

2 = rotational frequency (cycles per sec.) 


Accordingly, 


bn = (wr? Win” fae" 


The values of b, obtained from the Myklestad re- 
sults were calculated and compared with the values 
given by Horvay. Frequency values required in the 
calculations were read from Fig. 5, using the plotted 
stationary values and the value at 194 r.p.m. The 
results obtained are given in Table 3, together with the 
comparable values given by Horvay. 


TABLE 3 


Myklestad 


No. 1 blade Vo. 2 blade Vo. 4 blade Horvay 

by b 5.15 bh = 4.98 db; 5.62 
bs b 12.58 be = 11.65 be 12.6 
bs 22.63 b 23.55 bs 22.27 b; 22.4 
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PREDICTED RESONANT CONDITIONS 


The conditions under which resonant bending may be 
anticipated can be predicted from a resonance diagram 
of the type commonly used in vibration work.'* Fig. 5, 
which is typical of such diagrams, predicts resonant 
conditions for the three sets of Cornell Aeronautical 
Laboratory fiber-glass blades. Reference will subse- 
quently be made to this diagram, with the objective of 
showing that the resonant conditions predicted were in 
fact consistently and accurately obtained. 

Construction and interpretation of the resonance dia- 
gram are both simple. First, curves are constructed 
showing the variation of blade natural frequency, 
as ordinate, against rotational speed, as abscissa. 
Curves are plotted for each mode being investigated. 
Second, a family of straight lines is constructed, radiat- 
ing from the origin, with slopes equal to the harmonic 
order numbers. The intersection of any natural fre- 
quency curve with any order line locates a resonant 
point. The order of the resonance is equal to the slope 
of the order line. 

Study of the resonance diagram leads to several inter- 
esting and important generalizations, which will be 
used later to correlate the results obtained from the 
analysis of the measured resonance data: 

(1) The frequencies at which resonance is predicted 
in any mode separate into groups. The groups do not 
overlap and are independent of blade stiffness. For 
example, it is seen that at all operating rotor speeds the 
first-mode frequencies fall between the second and third 
orders, the second-mode frequencies between the fourth 
and fifth orders, and the third-mode frequencies be- 
tween the sixth and ninth orders. 

(2) For even the exceptionally wide range of stiffness 
covered by the Cornell Aeronautical Laboratory blades, 
the variation of natural frequency with stiffness is small 
compared to the variation with mode. 

(3) A decrease in blade stiffness in any mode will 
always tend to decrease the resonance order, the effect 
of the blade stiffness becoming increasingly important 
in the higher modes. 

(4) An increase in rotor speed will, for all modes, 
always tend to decrease the resonance order, the effect 
of rotational speed becoming increasingly important in 
the higher modes. 

Referring to Fig. 4, it is interesting to determine how 
well the eighth-order bending oscillation shown can be 
correlated with the resonance diagram. Entering the 
diagram at the ground run-up speed recorded (167 
r.p.m.), it is found that eighth-order resonant bending 
of the medium-stiffness blade is exactly predicted. It 
is probable that the exact agreement obtained is some- 
what a matter of chance, since it is improbable that the 
computed natural frequencies are sufficiently accurate 
to provide such precise correlation. Nevertheless, the 
check obtained is extremely convincing. 
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Fic. 7. Sample oscillogram, similar to Fig. 6, for the medium-stiffness blade 
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Fic. 8. Sample oscillogram, similar to Fig. 6, for the high-stiffness blade 



































HARMONIC ANALYSIS OF BLADE BENDING MOMENTS reduced to the form of Fourier coefficients and phase 
: , : ‘ angles. 

In the course of the variable-stiffness program, careful 
consideration was given to the problem of evolving ; 
ne ; ; ‘ Data-Analysis Procedures 
methods of analyzing and evaluating data which would 
be convenient and accurate and, at the same time, The amplitude coefficients and phase angles were ob- 
provide an insight into the data not obtainable directly — tained from harmonic analysis of selected cycles of each 
blade beam bending trace for each test point. Each 


from the oscillograph records. 
trace was reduced to the first ten terms of its Fourier 


Since the rotor is a rotating mechanical device and 
the forces applied to the rotor in steady-state flight are expansion, using a Dent-Draper Rolling Sphere Ana- 
periodic, the possibility of reducing the measured data _ lyzer, of the Corradi type. The harmonic analyzer per- 
tothe form of a Fourier expansion appeared advanta- forms integrations that measure the sine and cosine 
geous. Reduced in this manner, the data could be components of the Fourier series expansion of the curve 
readily presented in tabular form. More important, being analyzed. From the components, the resultant 
however, was the fact that the reduced and tabulated amplitudes and phase angles are computed. The ex- 
harmonic data could be directly interpreted in terms of — pansion was carried through ten terms since preliminary 
such important parameters as number of blades, aero- analyses of shakedown data had shown that bending 
dynamic loads, blade bending modes, and similar effects in the second and third bending modes would pro- 
physically definable quantities. For example, the num- duce harmonic components, respectively, Letween the 
ber of blades determines the phase relation between the fourth and fifth and sixth and ninth orders. 
loads applied to the helicopter; aerodynamic loads are Conversion of the harmonic coefficients to units of 
calculated in the form of summations of harmonically blade beam bending moment was accomplished by 
varying quantities having frequencies that are multiples applying a channel sensitivity factor, determined from 
of the rotor speed, and, as revealed in this paper, cer- dead-weight calibrations recorded before and after 
tain important bending characteristics of rotor blades each flight. These values were further adjusted for 
can be readily explained in terms of harmonic bending — sensitivity and drift variations, which in some cases 
components. For these reasons, all rotor data were occurred during the flight, by comparing calibration 
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Fourier reduction of transition bending data (Fig. 6) for the low-stiffness blades, showing predominance of fourth-order 


bending in the second mode and sixth-order bending in the third mode. 


injects taken at intervals during the ground calibration 
with injects taken between the data record runs. Calli- 
bration injects were obtained for each trace by sequen- 
tially switching the data strain-gage bridge for each 
channel to a calibrated signal of fixed value. 

Sample oscillograms showing sections of the flight- 
test records obtained for three test points are repro- 
duced in Figs. 6, 7, and 8. The nine channels recording 
beam bending data may be identified from Fig. 3. 
Three traces on each record have been “‘heavied up”’ so 
that they may be followed. The three traces (Station 
72 shown as a dash line, Station 135 as a dot-dash line, 
and Station 234 as a solid line) were chosen to bring out 
certain facts concerning the bending characteristics of 
the blades. The manner in which they do so can be 
visualized by reference to the mode shape diagrams 
(Fig. 3). The diagrams show that Station 72 is located 
so as to exhibit bending content in all three modes. 
Station 135 is located near a point of inflection in the 
second mode and should therefore show little second- 


mode content but large first- and third-mode content. 
Station 234 should show little first-mode content be- 
cause of its outboard position and substantial second- 
and third-mode content because of its position between 
the second- and third-mode inflection points. 


The sample records, all taken in the transition range, 
exhibit clearly the predominant higher mode bending 
components. Visual inspection of the records for the 
high- and medium-stiffness blades (No. 4 and No. 2 
reveals, respectively, large fifth- and fourth-order com- 
ponents. Similar inspection of the record for the low- 
stiffness (No. 1) blade reveals an even larger fourth- 
order component. It is seen that the harmonic order 
and spanwise distribution of the higher frequency com- 
ponents agree with the concept of higher mode bending. 
Specifically, it is seen that the second-mode bending 
component, which appears as a fifth-order component 


with the No. 4 blades, appears as a fourth-order com- 
Further, with the No. | 
blades, which provide a further decrease in stiffness, the 


ponent with the No. 2 blades. 
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fourth-order component becomes increasingly large, 
apparently denoting a closer approach to resonance. 
By comparing the phase of the bending traces for Sta- 
tion 72 and Station 234, it is seen that the 1S0-deg. 
phase difference required between bending measured at 
blade stations on opposite sides of the second-mode 
inflection point exists. 

While the particular records shown clearly exhibit the 
characteristics, oscillatory amplitudes, and harmonic 
relationships associated with resonance, it is to be 
noted that they show extreme conditions and that the 
records for many other test conditions reveal, to the 
eye, little or none of the pronounced resonance visible 
here. On the other hand, as will be shown later, the 
harmonic analysis of almost every record does disclose 


substantial higher harmonic content. 


Fourier Moment and Phase-Angle Plots 


The harmonic data derived from each of the preced- 
ing sample oscillograms are shown in Figs. 9, 10, and 11. 
The ten amplitude coefficients and phase angles for 


Fourier reduction, analogous to 
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Fig. 9, for the medium-stiffness blades 


each of the nine bending data channels are plotted 
against blade span. For clarity, harmonics one 
through five are shown on one plot, and harmonics six 
through ten on another. It is interesting to examine 
the harmonic plots and to determine how they may be 
identified with the oscillograms and interpreted in 
terms of the resonance diagram. 

From a resonance point of view, the harmonic data 
for the high-stiffness (No. 4) blade reveal several impor- 
tant facts. It is seen (Fig. 11) that the predominant 
harmonics are the first, third, fifth, and eighth. The 
strength of the first harmonic is explained on the basis 
of the usually strong first-harmonic air loading. The 
fifth and eighth harmonics, which cannot be explained 
on the basis of known air loads, are shown by the reso- 
nance diagram, Fig. 5, to fall on predicted resonant 
points. Their large response seems thus explained. 
It is also seen from the resonance diagram that the 
large amplitude of the third harmonic appears to be 
explainable on the basis of resonant amplification, since 
the high-stiffness blade approaches resonance with the 
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third-harmonic loads. Referring to the oscillograph 
traces, Fig. 8, it is seen that the large fifth-order com- 
ponent is readily identified in the trace (No. 1) for Sta- 
tion 234. 

The spanwise distribution of the harmonic amplitudes 
and phase angles is next examined in order to determine 
whether they are consistent with the blade bending 
mode shapes. Referring to the third-harmonic ampli- 
tude and phase plots of Fig. 11, it is seen that the span- 
wise distribution of moment approximates the first 
bending mode shape in that it decreases smoothly from 
a peak value at mid-span to zero at the blade tip and 
that the peak amplitudes are all reached at approxi- 
mately the same azimuth position. Referring to the 
fifth-order amplitude, it is found that the spanwise 
distribution decreases from a peak value at 0.25 span, 
to nearly zero at 0.56 span, then builds up to a second 
peak at 0.75 span, and drops off to zero at the tip. 
From the phase-angle plot, it is seen that the peak 
values for the inner half of the blade reach a maximum 
at approximately the same azimuth position (15°) and 
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to Fig. 9, for the high-stiffness blades. 


that the peak values on the outer half of the blade are 
all reached a half cycle later (180°). It is seen, there- 
fore, that the second-mode amplitude distribution is 
consistent with the second-mode bending shape and 
that second-mode bending is verified. Similarly, by 
noting the two cusps in the eighth-order bending ampli- 
tude distribution and the double phase shift from 
approximately 30 to 210 to 30 deg., large bending in the 
third mode may be verified. 

Study of the harmonic data for the medium- and 
low-stiffness blades will similarly establish the existence 
of large bending contributions in the second and third 
modes and show that the occurrence and order of the 
resonances measured are consistently predicted by the 
resonance diagram. 

It is particularly interesting to note the extremely 
large fourth-order component apparent in the sample 
record, Fig. 6, for the low-stiffness blades. As _ seen 
from the record and shown in the harmonic plots, Fig. 9, 
the contribution of higher mode bending (fourth and 


sixth orders) is larger than the sum of all other har- 
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Fic. 12. Envelopes showing spanwise distribution of maximum 
revolution. 
monics. Although no substantiating evidence is given 


here, it is noted that the peak bending moments re- 
corded for this test point are, over most of the blade, 
more than one and one-half times as great as the de- 
sign moment envelope obtained from the maximum 
moments computed for any and all flight conditions. 


ADVANCE RATIO EFFECTS 


The data obtained from the test program show that 
the critical blade bending conditions were not experi 
ences in high-speed flight but rather in the transition 
region. Further, it was generally found that a sub- 
stantial, and sometimes predominant, portion of the 
total bending moment consisted of third- and higher 


order components. 


This section of the paper presents results that sub- 
stantiate the foregoing statements. Two types of data 
plots are reviewed. The first type (Fig. 12) consists of 
bending moment envelopes for the three transition 
conditions shown by the sample oscillograph records, 
plotted against blade span. The second (Figs. 13, 14, 
and 15) consists of bending moment envelopes plotted 


against advance ratio. In the latter plots, data for the 


and minimum beam bending moments measured during one rotor 
The envelopes correspond to the oscillograms shown in Figs. 6, 7, and 8. 


nine blade stations are condensed on a single figure, 
thus providing a convenient means of studying the 
spanwise variation of the advance ratio effect. 

Current methods of blade stress analysis assume that 
the total moment applied to the blades consists of a 
combination of only steady, first-harmonic and second- 
harmonic bending moment components. In order to 
provide data that might be used to study the reliability 
of these methods, one set of envelopes including only 
the first- second-harmonic the 
measured data was prepared. This may be compared 
with the theoretically predicted harmonic components 
A second set 


and components of 


for the corresponding flight conditions. 
of envelopes has also been determined for the sum of 
ten harmonics. The two envelopes, plotted on the 
same figure, provide a convenient method of showing 
the measured contribution of the higher harmonic terms 
and of disclosing how large they become in certain 
conditions. 

By calling attention to these results, it is hoped that 
the structural importance of the resonance phenomenon 
throughout the speed range will be established and that 
attention will be called to the need for better design 


criteria. 
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Spanwise beam bending moment distributions for the low-stiffness blade, showing the contributions of the sum of third 


and higher order bending terms, as a function of advance ratio. 


Bending Envelope Synthesis 


The moment envelopes presented are plotted from 
data computed with a ‘‘synthesizer’’ developed for use 
in the program. The synthesizer is an electromechani- 
cal computer that provides the sum of either the first 
two or first ten harmonics, derived from the measured 
data, as a function of blade azimuth. The data set 
into the machine consist of the Fourier harmonic co- 
efficients and phase angles for each trace. 

While normaily the envelopes might be determined 
from peak moment values measured directly from the 
oscillograph records, this could not be done in the 
present case without sacrificing accuracy because filters 
inserted in the data channel circuits for the purpose of 
removing high-frequency hash also produced substan- 
With the filters used the re- 


? 


sponse at frequencies as low as 3 cycles per sec. (rotor 


tial signal attenuation. 


speed) fell off to 90 per cent, and, at frequencies corre- 
sponding to the important fourth and fifth harmonics 
(12 to 15 eycles per sec.), the response fell to 75 per 
cent. The filter attenuation problem was one of the 
reasons for employing the synthesizer; another was 
that use of the synthesizer also provided a convenient 
means of obtaining the harmonic sums used to construct 
the moment envelopes for the sum of the first and 


second harmonics 


The steady bending moment distribution about 
which the loading oscillates throughout the azimuth 
cycle is also shown in the maximum-minimum envelope 
The steady moments were computed from the 


Adding the maxi- 


plots. 
test results, using the synthesizer. 
mum and minimum cyclic moments to the steady 
moment determined for each blade station gave the 
maximum and minimum blade beam bending moment 
envelopes shown in the figures. 


Description of the Harmonic Synthesizer 


Basically, the harmonic synthesizer is a device in 
which the coefficients, phase angles, and harmonic re- 
lationships of a Fourier series may be inserted using a 
combination of electrical and mechanical means and 
the values of the sum of the series determined as a func- 
tion of angular displacement. The machine consists 
principally, of ten pairs of selsyn units, ten potentiom- 
eters, and ten sets of gears having the exact ratios of the 
integers one through ten. The proper phase relation- 
ships of the harmonic inputs are obtained by setting 
the rotor positions of ten of the selsyns. The rotors of 
the remaining ten selsyns are geared to a common shaft 
by the ten sets of gears and are, therefore, caused to 
rotate through angles proportional to their harmonic 


order. The potentiometers provide a means of regulat- 
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ing the voltage-—proportional to harmonic amplitude 
applied to each pair of selsyns. By adding the voltage 
outputs of the ten pair of selsyns in series, as a func- 
tion of angular position of the common geared shaft, 
the vector sum of the ten inputs is determined. When 
applied to the present problem, this is equivalent to 
providing the sum of the harmonic bending moment 
components as a function of blade azimuth. The sum 
of the first and second harmonics was obtained by a 
switching arrangement that added the voltage output 
of only the two proper sets of selsyns. 

It was stated previously that the test data showed 
that the critical blade bending conditions were measured 
in the transition region. Also, it was shown by one 
example for each blade that a large part of the total 
moment consisted of the higher harmonic terms and 
that these terms originated from resonant bending. 
Since data are only given for one test point with each 
blade, the only comparisons that can be made directly 
are those concerned with blade stiffness effects. From 
the envelope data given in this section, the role played 


by advance ratio effects may be evaluated separately 


from the resonance effects. 


In Fig. 12, envelopes synthesized from the harmonic 
data given in the harmonic plots (Figs. 9-11) are 


Spanwise beam bending moment distributions, similar to Fig. 13, for the medium-stiffness blade. 


plotted against blade span. Each envelope shows the 
maximum and minimum moment experienced at each 
blade station during one rotor revolution. One set of 
curves, for the sum of ten harmonics, shows the total 
The second set of curves shows 
The steady 


moment at the station. 
the sum of the first and second harmonics. 
moment curve shown applies to the envelope for the 
sum of ten terms. Although the envelopes for the sum 
of two terms are constructed by plotting the maximum 
and minimum moments relative to the steady value, it 
should be recognized that the steady value does not 
have the physical significance with respect to two terms 
which it has with respect to ten terms. It is also noted 
that the steady value is not the average of the maximum 
and minimum values. This follows from the fact that 
the phase differences between the various harmonics 
cause the maximums and minimums to occur at dif- 
ferent azimuth positions and to sum to different total 
values. In certain cases, it is found that the effects of 
the phase differences will even cause the sum of only 
the first and second harmonics to be greater than the 
sum of the ten harmonics. 

The important bending moment contribution of the 
higher harmonics in the transition region is apparent in 
Fig. 12. It is seen that, in addition to accounting for a 
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large part of the total moment on all blades. the higher 
harmonics may produce substantial negative bending. 
This is found to be particularly so with the low-stiffness 
blades, the negative bending resulting largely from the 
predominant fourth-order second-mode bending com- 
ponent. 


The effect of advance ratio can be shown by plotting 
the maximum-minimum moment against advance 
ratio, as has been done in Figs. 13, 14, and 15. AI- 
though curves are shown for only the normal rotor 
speed, forward center of gravity, level flight conditions, 
results for other conditions of operation show that 
these curves are typical. In the form shown, two 
main effects are observed. First, at all advance ratios 
and all blade stations, the higher harmonic contribu- 
tions are important. In most cases, and particularly 
with the medium-stiffness blades, the envelopes peak 
twice, first in transition and later at high speed. In all 
Over the in- 
board half of the blade the contributions of the higher 
harmonics are large in the transition conditions. In 


cases the transition peak is greatest. 


high-speed flight, the higher harmonic contributions 
continuously decrease although the peak values are 
again increasing. The small differences between the 
sum of two and sum of ten terms obtained over the 


Spanwise beam bending moment distributions, similar to Fig. 13, for the high-stiffness blade 


inboard half of the blade in high-speed flight do not re- 
sult from absence of the higher harmonics but from the 
fact that the predominant harmonics are out of phase 
and cancel one another. Over the outboard half of 
the blade, the contribution of the higher harmonics at 
intermediate speeds tends to be larger than inboard. 
This can again be accounted for primarily by a phase 
effect, since the large second-mode contribution has 
Also, in this 
region of the blade, the second- and third-mode con- 


shifted phase by approximately 180°. 
tributions are in phase. At high speeds, the differences 
between the harmonic content curves become extremely 
small, and, in some cases, the curves cross. This re- 
sults from the fact that near the blade tip the second- | 
and third-mode components are out of phase and that 
the amplitudes of the third-mode components become 
larger than those for the second mode. 
‘ 
RESONANCE CORRELATION | 
In the case of the three transition test points pre- 
viously analyzed, it was shown that the occurrence and 
order of the higher harmonic bending components 
agreed well with those predicted by the resonance dia- 


gram. It is shown in this section that equally good 


agreement can be obtained by treating the 97 test 








point: 
whole 
of eat 
speed 
the e! 
and f 
are a 
condi 
all di 
194 1 
seen 
with 
point 
Th 
have 
appli 
bar | 
coml 
rotor 
age I 
were 
of it 
and 
pape 
It 
poin 
reso! 
as fe 
mod 
seve 
ings 
can 
the 
li 
reso 
tion 
agre 
wer 
prec 
(Fig 
Tha 
cree 
tha 
shif 
first 
tior 
tior 
fiftl 

















RESONANT BENDING AND 
points recorded in the test program (Table 1) as a 
whole. This is done by determining the average value 
of each harmonic for all test points in which the rotor 
speed and blade stiffness are the same. In this way 
the effects of advance ratio, center of gravity position, 
and flight condition (level flight, climb, autorotation) 
are averaged out of the results. As an example, 19 test 
conditions were averaged to obtain the average value of 
all data for the medium-stiffness blades operating at 
194 r.p.m. Referring to Table 1, the 19 points are 
seen to consist of seven points with forward and five 
with aft helicopter center of gravity positions and seven 
points in climb. 

The averaged values obtained for three blade stations 
have been plotted in Figs. 16, 17, and 18. Each figure 
applies to a particular blade spanwise station. Nine 
bar graphs are plotted on each figure, one for each 
combination of three blade stiffnesses and three test 
rotor speeds. The three blade stations for which aver- 
age harmonic plots are given (Stations 72, 135, and 234) 
were selected because they provide a maximum amount 
of information with which to correlate the predicted 
and measured resonances, as discussed elsewhere in this 
paper. 

It is recalled that in an earlier discussion 
pointed out that the harmonic orders predicted for 
resonance in each mode fell into nonoverlapping groups, 
as follows: first mode—second or third order; second 
mode—fourth or fifth order; 
seventh, eighth, or ninth order. 
ings in mind, the predominant harmonics in each group 
can be selected from each bar graph and tabulated with 


and third mode—sixth, 


Keeping these group- 


the results shown in Table 4. 

Inspection of the table clearly reveals that rotor 
resonant bending is a phenomenon common to all condi- 
tions of helicopter flight and that the measured data 
agree consistently and accurately with the results that 
were predicted. It is noted that each and every trend 
predicted from inspection of the resonance diagram 
(Fig. 5) is systematically and accurately verified. 
That is, the orders of the harmonics decrease with de- 
crease of blade stiffness and increase of rotor speed and 
that the harmonics associated with the third mode 
shift order more often than those associated with the 
first or second mode. Finally, it is seen that the loca- 
tion of Station 135 near a point of second-mode inflec- 
tion is reflected in the absence of significant fourth- or 
fifth-order harmonic content and that the predominant 
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content at Station 72, averaged from representative data for all 
régimes of helicopter flight 


second- and third-mode components anticipated at 
Station 234 are obtained. 


CONCLUDING REMARKS 


It is obvious from the results presented that the blade 
resonance problem is of primary concern to rotor-blade 
designers and that it should be carefully considered in 
the design and flight evaluation phases of rotor develop- 
ment programs. 

While the good correlation obtained between pre- 
dicted and measured results indicates that it should be 
possible to predict resonant conditions during the de- 
sign stage, unfortunately, it is found that little can be 
done to avoid coming close to one or another resonant 
point. It appears that the best that can be done, at 
least until the nature of the resonance is understood, is 
to make provision in the design so that the natural fre- 
quencies do not coincide with multiples of the normal 
rotor operating speed and, more important, to provide 


Station 72 


Tabulation of Predominant Harmonic Orders 


TABLE 4 


-Station 135 


Station 234 


175 194 P06 175 194 206 175 194 206 
No. 4 2: 5:9 3, 5,8 3. 5,8 3 ,9 3, ,8 3, ,8 SS. 1,8 » 4,8 
No. 2 a ae a. 4.7 a. 4,7 2 ,3 2,—,7 S; 7 oe | »4,7 
No. 1 2, 4,7 2, 4,6 2. 4.6 2 + 2,—,6 :. ,6 4,7 » 4,6 » 4,6 


Dashes indic: 


ite 


> that the content is 


relatively 


small. 
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Fic. 17. Bar graphs showing harmonic bending moment content 
at Station 135, averaged from the same test points as Fig. 16. 
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Fic. 18. Bar graphs showing beamwise bending moment content 
at Station 234, averaged from the same test points as Fig. 16. 
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sufficient structural strength, over and above that re- 
quired to carry the normally computed steady, first. 
and second-harmonic loads, so that the blade can safely 
operate subject to higher order loads and resultant 
moments. <A tentative but detailed method for esti. 
mating higher order loadings, based on the results of 
the variable-stiffness program, has been outlined by 
Daughaday and Kline.'* 

Anticipating the probability of higher order bending 
stresses, it is important that the blade designer recognize 
the critical role that resonant conditions may play in 
determining the fatigue life of blades. Since the blade 
resonant frequencies are harmonic multiples of the rotor 
speed, it is important to anticipate that the number of 
load cycles applied to the blade will be proportional to 
the product of total rotor revolutions and the excited 
harmonic orders and that the peak stresses and vibra- 
tory moment amplitudes will be larger than predicted 
by current theory. 

Finally, the results call attention to the necessity of 
providing sufficient data stations, when conducting 
blade strain surveys, so that information will be ob- 
tained at all critical blade locations. Equally impor- 
tant, it is necessary that the stations be judiciously 
chosen so that they do not fall at points of inflection 
and thereby fail to detect important higher mode con- 
tributions. 

It is trusted that this paper will direct the attention 
of those concerned to the rotor-blade higher mode 
bending phenomenon and that it discloses the most 
important manifestations of the problem. It is recog- 
nized that no answers are provided to many of the ques- 
tions raised. Such answers are not yet available. It 
is hoped, however, that this paper will occasion suffi- 
cient interest and concern so that the problem will re- 
ceive the attention it merits and that the needed answers 
will be forthcoming soon. 
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Further Note on the Use of the Neutral Point as 
a Stability Parameter 


S. J. Deitchman 
Aerodynamicist, Bell Aircraft Corporation, Buffalo, N.Y. 
March 11, 1953 


I* A RECENT NOTE in the Reader’s Forum of the JOURNAL, Baron! 
pointed out that the neutral-point location is not necessarily 
au indication of longitudinal static stability if the control effec- 
tiveness changes with angle of attack due to downwash effects. 
In showing this, it was assumed that the control effectiveness is 
independent of control deflection. For the type of aircraft dis- 
cussed, however (canard, supersonic missile), this assumption is 
excessively restrictive. 

The nonlinearity of the downwash at the aft wing derives from 
the behavior of the vortex sheet behind the forward wing. If the 
forward wing span is large with respect to the aft wing and the 
tail length is small with respect to the forward wing span, the 
trailing vortex sheet can be considered essentially flat and the 
downwash variation with normal force linear. However, with 
the forward wing of approximately the same span and a relatively 
large distance ahead of the aft wing, the trailing vortex sheet is 
partially or completely rolled up at the aft wing location. The 
vortex cores from the tips will either cross the aft wing or pass 
close to its tips. Therefore, the distribution of downwash will 
be nonuniform, and, furthermore, it will vary with Cy since the 
degree of rolling up of the trailing vortex sheet at the aft wing loca 
tion depends on Cy.” Since deflection of the elevator on the 
forward wing represents a change of Cy of that wing, the down- 
wash effect on control effectiveness may be expected to depend 
ou control deflection as well as a. It should be noted that the 
use of the trailing vortex concept for supersonic downwash has 
been justified by various authors. 4 

Following the pattern and using the notation of reference 1, but 
using the NACA convention for writing derivatives, the sta 
bility equation may be derived in a more general way. Let the 
variation of pitching moment and normal force coefficients with 
control-surface deflection at constant a be of the general nature 
shown in Fig. 1. 

Then Cn, 
efficient due to the variation of Cg and Cy; with 6, and the equa- 


and Cy, will be increments of moment and lift co- 


tion for the pitching moment about a point } in terms of the con- 


ditions at point ais 


(Culp = (Cu) — CyuX* + [(Cn,), — (Cn,)X*)6 + 
oO 0 


fn(6) — fr(6)X* (1) 


where Cy. = fm(6) and Cy, = f,(6) are defined in Fig. 1. 


mo 

From this, differentiating with respect to a, determining (Cm a 
in terms of conditions at the neutral point, .V, letting Xa,,* go to 
zero at the center of gravity and Yy* be the distance from the 


Contributions should not exceed 800 words in length 


+ 


c.g. to the neutral point (positive forward), and dividing through 
by Cya, the stability equation becomes, for trimmed flight, 


(Cagyess? = Xn* [ + dy in. (cx) | + 
- : d ; d 
(Sr — 5x) 7 (mse © oe SES Oe 2 
where 
F(é7) = fnldr) — fr(b7)X* = fn(br) 
F(6v) = fnl6n) — fildn)X n* 


For the missile to be stable, the right-hand side of Eq. (2) must 
be negative. Whether or not this will be the case depends on the 


magnitudes and directions of the variations of ( ns and Cy,, as 
é 5 


well as the relative locations of the c.g. and neutral point 
For'a specific example, if the pitching-moment characteristics 
» 


of the missile are as shown in Fig. 2, then the following relations 
will probably hold for trimmed flight : 


bv <0 
(67 — by) Q if Yn" 0 
d ( d . 
v;) <0 — (Cmz) 0 
dCy “ dCy lac 
from which 
ds. 
(67 — by) - (Ce leg 0 
dCy : 
if XN U 
d 


dCs [F(ér) — F(éy)| 2 0 | 


and 
bxn(d/dCn) (Cy;) = 0 


where the upper signs or lower signs of 2 in paired conditions are 
to be used simultaneously. 
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Applying these inequalities to the right-hand side of Eq. (2 
shows that, if certain conditions of the magnitude and direction of 
pitching moment variations with 6 are satisfied, it is possible to 
have a stable missile with the c.g. aft of the neutral point or an un 
stable oue with the c.g. forward of the neutral point. Thus the 
concept of static margin as an all-inclusive measure of stability 
breaks down. 

The chief significance of this result is the difficulty it intro 
duces into making preliminary estimations of missile stability 
If experimental data are available, the degree of stability is still 
indicated by the slope (Cen Je.g.. 7, but the instantaneous aero- 
dynamic center location of the complete configuration is no 


longer given by this slope. 
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A Note on the Weight of the Blades in an Axial- 
Flow Compressor 


John Mitchell Stephenson 

Research Associate, Division of Engineering, Brown University, 
Providence, R.I. 

March 19, 1953 


¢ bs GENERAL PROBLEM of designing an axial-flow compressor 
for the minimum weight is complex and can usually be 
discussed only for particular numerical cases. But by making 
some simplifying assumptions—especially that all the stages in a 
compressor are aerodynamically similar—it is possible in a few 
minutes to calculate the weight of all the blades in terms of the 
weight of the blades in the first stage. 

A chart can then be drawn up to show the weight penalty in- 


curred if manufacture is simplified by using the same blades 
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(suitably shortened) for more than one row. The only proper- 
ties of the compressor that must be known are the number of 
stages and their pressure ratio and efficiency. It is assumed 
that the maximum allowable bending stress is independent of the 
centrifugal stress. If this is not admitted and only the maximum 
total stress is fixed, the weight penalties for grouping are slightly 
larger than those given here. 

The weight W of a row of blades is proportional to their num- 
1, their height h, and a factor f, 


Thus 


ber B, their root section area / 
which takes account of the taper. 


Wa BAhf 


But, approximately, A « ct, where cand?¢ are the chord and maxi 
mum thickness of the root section and where the product Bc is 
proportional only to the radius and solidity of the root section 
Therefore, if we assume that Bc, f, are the same for all stages, 


W« th (1) 


It was shown in reference (1) that the maximum bending stress 
in a blade can be written approximately as o, where 
ot?/h = (cos ¢ + Bsin ¢)-QC,AT/2uBck (2) 
where (cos ¢ + 8 sin ¢)/k is a number of the order 10, depending 
of the geometry of the blade; Q is the mass flow; C,AT is the 
work done in the stage; and uw is the mean tangential blade 
speed. The right-hand side of Eq. (2) depends on the aerody 
namic design of the stage, and we shall assume it to be the same 
for all of the stages in a compressor. Hence, if the stress is also 


constant, 


and substituting in Eq. (1) 
W« h'?? (3) 
Let us also assume that the Mach Number is the same at the 
entrance to each stage. Then Q»y/7/PA is constant through 
the compressor, where 7 and P are the temperature and pressure 
(either total or static values will do) and A is the flow area (annu- 
lus area X cosine of the flow angle). Therefore, if the mean 
radius is constant, 
h<« A«+/T/P 
and, from Eq. (3) 
Wea T*/s. p-’ (4) 
Now assume that the polytropic efficiency » is constant through 
the compressor and that the entry temperature and pressure are 
fixed. Therefore, 
Pa« TYV(% 1 


and, if y = 7/5, 
P« T? 5n 
Any appropriate value for 7 can be used; a convenient value is 
n = ®/;, which gives 
Hence, Eq. (4) becomes 
Wa P- (5) 
If the pressure ratio of each stage is the same, the pressure at 
entry to the mth stage is 
P,« r*-! 
and the weight of the mth stage is, from Eq. (5), 
y,« £,-~ 
or 
W,« exp [—(5/4) (nm — 1)-log, r] (6) 


So the stages can be made progressively lighter. 
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The total weight of all the blades is clearly least if those in 
each successive row are thinner as well as shorter, so as to keep 
the stress constant at the maximum allowable value. But it is 
frequently cheaper to manufacture successive rows of rotors (or 
stators) by cropping the tips of the rotor (or stator) blades in the 
previous stage, thus keeping the root dimensions constant for 
two or more stages. In this case, the stresses in the second, third, 
etc., rows in each group are less than the maximum allowable 
value, and there is a weight penalty. 

If there are ¢), go, gs, .. ., blade rows in each group and N rows 
altogether, the total weight is 2 times the weight of the first stage 


row, where 


or, writing [(5/4)log, r| as g in Eq. (6), 
2 = gi + go exp (—qgi) + gs exp [—q (gi + g2)] (7) 


‘he extreme conditions are when the stages are all the same and 
Q = N, and when they are all different, so that 


@ = [1 — exp (—qN)]/[1 — exp (—q)] 


Between these limits, the condition that 2 is a minimum for a 


given number of groups is 


02 on ov? 


= ~ a = () 
Og: Og» = O8's 
or 
1 — qg2 = exp (qg:)) 
1 — qgs = exp (qg2)} (8) 
1 — gg, = exp (qgs) 
ete. 
The application of Eqs. (7) and (8) can best be illustrated by a 
numerical example. Suppose g = '/s; then using Eqs. (8) we 


can plot, as in Fig. (1), curves of go, (gi: + ge), gs, (g1 + g2 + gs), 
£4, (£1 + go + gs + gs), etc. Then if we have a compressor with, 
say, eleven stages, we draw vertical lines intersecting the lines of 
(g1 + go), (gi + ge + gs), and (g: + go + gs + gy), at the ordinate 
eleven. The appropriate values of g), ge, g3, and gy are read off 
the intersections to the nearest whole number. Therefore Table 
1 can be compiled, the corresponding values of 2 being calcu- 
lated from Eq. (7). The weight penalty incurred by grouping 
is seen at a glance in the last row. The relative weight of dis 
tributions other than those given in the table can easily be cal- 
culated from Eq. (7). The case of two groups, for instance, is 


shown in Fig. 2. 
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TABLE | 
Number of Groups ] 2 3 4 11 
Distribution 11 4-7 3-3-5 2-2-3-4 
Q2 6.58 5.53 5.02 1.2] 
Q 31 1.19 1.00 


4.21 2.61 1. 56 I 


In practice, the apparent advantage of using many groups may 
not always be realized because of a lower limit to the blade thick 
ness imposed by the process of manufacture. 
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A Method for Differentiation of Experimental 
Data 


William M. Kaufman and Marvin Shinbrot 
Ames Aeronautical Laboratory, NACA, Moffett Field, Calif 
February 13, 1953 


i ANALYZING THE MEASURED TRANSIENT RESPONSE of aircraft 
to disturbances, it is often desired to evaluate time deriva 


tives of the measured quantities. Among other uses, these 


quantities can be employed in various forms of the derivative 


method for determining the aircraft stability parameters.! How 


ever, there has been a tendency to avoid this procedure because 


of the difficulty in evaluating the time derivatives with sufficient 
accuracy. 
Receiutly, one of the authors was faced with this problem while 


attempting to develop methods for analyzing the measured 
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transient response of nonlinear systems. The usual graphical 
and numerical differentiation methods were checked first by 
applying them to analytically expressed functions having certain 
characteristics of available or expected test data. These meth 
ods did not consistently yield even the first derivative to the de 
sired accuracy when the basic data were rounded off to simulate 
recorded data typical of dynamic flight tests. This difficulty led 
to the development of the following differentiation method, 
which has been applied successfully to several problems at hand 

Basically, it is a differencing method refined by a curve-fairing 


procedure. Given a curve x = f(t), consider any point ¢ = T 
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Calculate 
Ax, = x(T + At) — x(7 
Ax. = x(T + 2At) — x(T) 
Ax-, = x«(T — At) — x(T) 
Ax-2 = x(T — 2At) — x(T) 


where Af is a convenient small interval. Next find the quanti 
ties Ax, /RAt(k = 1, 2, —1, —2), which are all approximations, 
of the familiar general form Ax/At, to the derivative dx/dt at 
t = 7. Then plot Ax,/kAt as a function of k, and fair a curve 
through the four points. The ordinate of this curve at k = 0 is 


a measure of dx/dt, since 


x(T + kRAt) — x(T) , x(T + At) — x(T 
lim = lim 
k—o RAt At—>oO Al 
dx /dt 


More than four points can, of course, be used but generally have 
not been found necessary; two points yield the often crude 
chord-line approximation 

This procedure is illustrated in the following test problem 
The function 


x = 6sin¢ + 0.5 sin 3¢ 


(0 < ¢t < wm) was selected as a troublesome one having a term 
(0.5 sin 3¢) with a relatively small effect on x but successively 
larger effects on dx/dt and d*x/di*. This function is plotted in 
Fig. 1. The only data furnished a computer who evaluated 
dx/dt and d*x/dt? by the above method were the values of x at 
Af intervals of 0.2 


or 8 in the second decimal place, corresponding to recorded data 


These x values were rounded off to 0, 3, 5, 


with a readability of about 0.5 per cent of the maximum value 
Plots of Ax; /kAt for a few values of 7 are given in Fig. 2. The 
resulting values of dx/dt (Ax/kAt at k = (0) are plotted as circles 
in Fig. 3. The entire procedure was then repeated with these 
dx/dt results to obtain the d*x/dt® points plotted as triangles in 
Fig. 3. Curves representing the true values of dx/dt and d*x/dt*, 
determined analytically from the given equation for x, are shown 
for comparison. The accuracy of the proposed method appears 
particularly noteworthy in view of the complicated shapes of the 
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derivative curves. 
d*x/dt? are 0.05 and 0.18, respectively—significantly smaller 
(especially for the second derivative) than those obtained when 
other well-known graphical and numerical differentiation meth- 


The root-mean-square errors in dx/dt and 


ods were applied to the given data. 
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Leading-Edge Separation from Delta Wings at 
Supersonic Speeds 


Mac C. Adams 
Cornell University, Ithaca, N.Y. 
March 23, 1953 


ie A RECENT PAPER, Legendre! has treated approximately the 
flat-plate delta 
His method, in brief, is to consider the problem as a po 


problem of leading-edge separation from 
wings. 
tential flow and to represent the separation by concentrated vor 
tices as in the figure. A further simplification is made by con 
sidering the flow in cross-planes perpendicular to the wing sur 
face (i.e., the Jones low aspect ratio theory). The restriction to 
supersonic speeds is inherent in Legendre’s treatment because 
conical flow is assumed. The method, however, may give an 
indication of the subsonic behavior if the aspect ratio is extremely 
small. 

In unpublished work at Langley Field (NACA) C. E. Brown 
and the writer have considered this problem using the same ap- 
However, there is a discrepancy in 
the determination of the lift, which will be briefly discussed be 


proach as that of Legendre. 
low. The purpose of this note is to give what is believed to be a 
correct interpretation of Legendre’s results and to show further 
that the simplified model of the real flow leads to some inconsist- 
encies 

Referring to Fig. 1, the more realistic model would have, instead 
of concentrated vortices, vortex spirals or sheets such as the lines 
A-B extending to vortex cones. Such lines would represent 
mathematical barriers across which there would be a potential 
jump but a continuous pressure. Now if the vortex sheets are 
not considered, as in the simplified model, the barriers A-B are 
Un- 


fortunately, however, the logarithmic singularity due to a con 


still necessary in order to render a single-valued solution. 


centrated vortex leads to the unrealistic result that a pressure 


jump exists in the field across the lines A-B. This same pressure 








jump exists across the wing surface in addition to the 
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pressure difference as computed by Legendre. By ignoring this 
pressure jump on the wing surface Legendre obtains extremely 
small negative lift coefficients for small angles of attack. By 
simply adding to Legendre’s result, the lift contributed by the 
pressure jump on the wing, as discussed above, the lift coefficient 
becomes, for small angles of attack, 


cs 1 
Cr = pa A.R.a oo Tv (: + <3) (4A.R y a 


where A.R. denotes aspect ratio and a the angle of attack in radi- 
ans. 

It is also interesting to calculate the lift by means of the well 
known impulse method of Munk. In this case a less lift results 
because the negative lift carried on the barriers A-B is also in- 
cluded. The lift coefficient given by this second method is, for 
small angles of attack, 


C, = (9/2) A.R.aw + x(4 A.R.)'/3 @’/3 


The inconsistency in these two different lift coefficients is clearly 
due to the simplified model considered (i.e., the condition of vor 
tex continuity is violated in ignoring the vortex feedlines to the 
cores ). 

It is not intended to say which, if either, of these results is cor 
rect. It is reasonable to say that the difference in the lift co 
efficients is a measure of the error introduced by the simplified 
model. It is noteworthy that both results give the same be- 
havior. The linear term is the same as Jones’ low aspect ratio 
result and the nonlinear term gives an increase in lift coefficient 
* The lift 
rectly for negative angles of attack. 


proportional to A.R a also changes sign cor 
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Improvement of Pressure Recovery of a Fixed 
Diffuser by Means of Suction 


Ashton M. Patterson + 

Institute of Aerophysics, University of Toronto, Toronto, Ontario, 
Canada 

February 13, 1953 


— RECOVERY in a supersonic wind tunnel is importait 
because the power requirements for a continuous tunnel de- 
crease and the running times of an intermittent wind tunnel in- 
crease with increase in pressure recovery. The pressure recovery of 
a fixed diffuser is limited by the ‘“‘starting condition”’ of the super- 
sonic wind tunnel. When the control valve of the intermittent 
blowdown wind tunnel is opened, the downstream pressure drops 
When 
) is reached in the throat of the nozzle 


and the flow of air through the tunnel is accelerated.! 
the velocity of sound (a+ 
(A«), the flow expands to supersonic velocity in the nozzle and 
is converted to subsonic velocity through a shock wave (starting 
shock wave). As the downstream pressure continues to decrease, 
the shock wave moves down the nozzle, and the velocity at the 
throat (Ax.) Once the 
(ax2) is reached in the diffuser throat, the diffuser acts as a nozzle, 


diffuser increases velocity of sound 


and no further drop in the downstream pressure will affect condi- 


tions upstream of the diffuser throat. The mass flow through 
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the two throats is 


pxdx Ax = px 212A «2 (1) 
° ° ° e . 
Using the energy equation and the relation a? = yp/p, Eq. (1 
can be changed to 
Ax/Ax2 = po,/ pi (2 


where (/o,/Po) is the stagnation pressure ratio across the normal 








shock wave. As the shock wave moves down the nozzle, it occurs 
it higher Mach Numbers, and the ratio po,/Po decreases in value 
and has a minimum value when the shock wave occurs in the test 
section (4,). As the starting shock wave will jump to some 
position downstream of the diffuser throat once it enters the dif 
fuser, the tunnel will start when 


Ax Arts S (po, Po) min (0) 


starting condition). The assumption that the wind tunnel 
starts with a normal shock wave introduces a small error into the 
calculation of the size of the diffuser throat that will allow the 
tunnel to run supersonically. [For short diffuser lengths the 
diffuser throat opening can be about 5 per cent smaller than that 
Therefore, Eq. (8) gives a good ap- 


“starting” 


> 


calculated from Eq. (3).] 
proximation for the calculation of the diffuser throat 
irea 

In order to improve the pressure recovery of the wind tunnel 
over that obtained with a diffuser having a throat area the size 
that the “starting condition”’ will permit, the diffuser throat 


If the throat area (Ax«-.) is smaller 


area must be made smaller 
than Ax /(po,/Po)min. then the shock wave will stand in the nozzle 
and the tunnel will not ‘‘start.”’ through the 


'Ax»and is equal to the mass flow through 


The mass flow 


diffuser throat is px2’ axe 
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the nozzle throat. (p«*2’ and ax.’ define the density and velocity 
of sound at the diffuser throat when the shock wave is standing 
in the nozzle.) If the shock wave was standing at the end of the 
nozzle blocks (starting condition), the mass flow through the 
diffuser throat would be pxod*2A*.. As the shock wave occurs at 


a higher Mach Number in the second case, px2 < px’. Therefore 
the diffuser throat cannot handle all the mass flow of the tunnel, 
and the mass flow that must be removed before the tunnel will 


run supersonically is 
m = (peo’ — pe2)dxeA «2 (4) 


If sufficient suction is applied upstream of the diffuser throat, 
the excess mass flow will be removed, and the tunnel will “start.” 
Once the tunnel is ‘“‘started,’’ the suction can be discontinued, and 
the wind tunnel will continue to run supersonically 

Preliminary experimental work was done in the 5- by 7-in 
blowdown-type supersonic wind tunnel at the Institute of Aero 
physics on the application of suction for initiating supersonic 
flow in the tunnel.2 The wedge-type diffuser used was con- 
structed so that slots of different shapes and sizes could be in 
serted upstream of the diffuser throat and so that the diffuser 
throat area could be adjusted between tunnel runs. The slots 
were connected by means of a flexible metal hose to a small 
vacuum tank, which was used as the suction source. The most 
efficient slot of those tested had a rounded inlet lip and a width 
the order of magnitude of the boundary-layer thickness.* With 
the diffuser throat area smaller than A*/(po,/Po)min-, the control 
valve was opened (the tunnel would not run supersonica!lly); suc 
tion was then applied, and once the tunnel started, the suction 
was discontinued 
[Eq. (4)] for the smallest 
19 cu.ft 
pressure of 127 mm. Hg and a density of 4.54 X 10 
had to be removed by suction before the tunnel would 

This amount checked well with the actual volume of air 


Calculation of the excess mass flow 
diffuser throat area obtained showed that per sec. at a 
‘ slugs per 
cu.ft 
start 
removed in the short time required for the tunnel to start (approxi- 
mately '/)9 sec 

Fig. 1 shows the improvement in pressure recovery that was 
obtained by using suction to initiate supersonic flow in the tun 
nel. The improvement in the pressure recovery was only of the 
order of 3 per cent over that obtained with a diffuser that would 
allow the tunnel to start without suction being applied. This 
small increase in pressure recovery indicated that the use of 
suction by itself would not give large increases of pressure re 
covery for fixed diffusers. Further work should be done on slots 
and suction systems to determine the maximum usefulness of 
suction It was found, however, that, if a long constant area 
throat section’ was inserted immediately downstream of diffuser 
throat (Fig. 2), the combination of the long throat and suction 
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to start the tunnel would improve the pressure recovery con- 
siderably. The improvement that was obtained in the 5- by 7-in. 
wind tunnel using a 6-in. constant-area throat length and suction 
was about 15 per cent over the best pressure recovery obtained 
with a fixed diffuser designed for the ‘‘starting condition.” 
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The Two-Dimensional Incompressible 
Aerodynamic Coefficients for Oscillatory 
Changes in Airfoil Camber 


Irvin N. Spielberg 

Wright Air Development Center, Wright-Patterson Air Force Base, 
Ohio 

March 2, 1953 


- gsr TYPES OF WINGS, notably those with low thickness 
ratio and low aspect ratio, are known to undergo rather pro- 
nounced deformations of chordwise sections when oscillating in 
their normal modes of vibration. Such deformations change the 
camber of the wing profile and, hence, can be expected to affect 
the aerodynamic pressure acting on the wing when it vibrates 
during flight. The significance, from the standpoint of flutter, 
of the nonstationary aerodynamic effects associated with camber 
changes is not fully known at present. In this note the two- 
dimensional, incompressible aerodynamic coefficients associated 
with harmonic changes in camber are derived. The derivation 
makes use of the linearized aerodynamic equations relating the 
pressure on the profile to the downwash distribution, as outlined 
in reference 1. Only one particular mode of camber deformation 
is studied—-namely, parabolic bending of the profile mean chord 
line—although the derivation can easily be extended to other 
chordwise deformation patterns. 

A thin airfoil of chord 2), immersed in a uniform stream of 
velocity V in the direction of the positive x axis, is considered to 
be undergoing small amplitude harmonic oscillations at right 
angles to the direction of the stream. Under the usual assump- 
tions of thin airfoil theory, the actual profile is replaced by its 
mean chord, which in the rest position is supposed to lie on the x 
axis extending from x = —1 (leading edge) tox = +1 (trailing 
edge) with x nondimensionalized by the semichord b. If z repre- 
sents the actual vertical motion of any point on the mean chord, 
positive down, then linearized theory gives for the downwash 


w = [(V/b) (02/0x)] + (02/dt) (1) 


Assuming the oscillations are simple harmonic of the form e'@', w 
being the circular frequency, Eq. 1 reduces to 


w = (V/b) [(0s/dx) + tks] (2) 
where k = bw/V is the reduced frequency. The subsequent de- 


velopment follows the procedure outlined in reference 1. A new 
chordwise coordinate, x = — cos @, is introduced and w is ex- 


panded in the Fourier series, 


wo, t) = Vel | Py +2 } 2 P,, cos no (3) 
1 


Then the pressure difference between the upper and lower sides 


of the mean line is given by 
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— 0 z sin #6 
T1(6, t) = — pV%e™ | 2a, cotan + 4 > a (4 
9 
2 1 n 


where, contrary to the notation of reference 1, the negative sign 
has been introduced here to give the pressure jump in the direction 
of positive forces—i.e., downward—in accordance with usual 
U.S. flutter notation. The following relations then hold between 
the coefficients P,, and a,,: 


a = C(k) (Po Pi) + P, j 
ik > > 5 
ay = (Pe—1 — Pet) af. 2x2 if 
where C(k) = F(R) + iG(k) is the Theodorsen circulation func- 
tion. 


Although this note is primarily concerned with the nonsta 
tionary aerodynamic coefficients associated with camber changes, 
with little additional work the coefficients for translational and 
pitching oscillations can also be obtained. It is therefore in 
tended to treat the following three types of motions: 


a(x) = he™, translation of the mean chord with 
amplitude h, positive down 

a(x) = ab[x + (1/2) ]et!, rotation of the mean chord 
about the quarter-chord point with angular | 
amplitude a, positive nose up p 

a(x) = ¢(1 — x?)e™!, parabolic change in camber of 


mean chord with amplitude ¢ measured at 
the mid-chord point, positive down 


The pressure distribution corresponding to each of these motions 
is obtained directly by applying Eqs. 2 through 5 as follows: 
Note that for convenience the factor ¢'“ is omitted temporarily 


For 21: 


w(@) = (iVk) (h/b) 


h 
Po =-(h)-, FP, =0. wn Zl 
b 
se Da h h 
a = [skC(k)]- = aor 
b b 
f kN h h (7) 
ay = ” = ai, 
( | b b 
dn =(), a>! 
For 2:: 


k 
w(@) = V (: +18) ik cos 6 | a 


2 Rk » k > 
Po = i oa ee a, Pi = t= Fie, P,=0, n>I1 
c k 
ao = |] (1 + 7k)C(k) — 1 5 | o& = da @ 
R2 
ay = (it — ) a = ig @ (8 
4 
ae = (k? tia = doq a 
An =, n>2 
For 2;: 
- V ; k on a] = k ~ 20 ¢ 
w3(0) = 1 9 + 2 cos i cs 2 ; 
Pp a _k ¢ ies ¢ 
Mo ; mee 
P, = ie =, P, =0 2 
4/ b 
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types 


where 


Substi 
tions J 


Ww, = 


bW3 = 


(a 


Ifn 
by b% 
defini 
expre: 


where 








’ 
(4 


ive sign 


lirection 
h usual 
set ween 


1 fune- 


1onsta- 
langes, 
al and 
re in- 


tions 
lows: 


ly 


( 
‘ 


é ¢ (9 
a = (ih) = aay 
y 

k7\ ¢ ¢ 
( = = ast 

8 J b ‘$ k 
a =(), mn 3 
Hence, 


y2}2 ’ ( h s 
-p 2 cotan doh t doa @ + Aor + 
2 b = 7% 


h y 
4 sin 6 (a + dia @ + ay + 2 sin 26 * 
b b 


( Pg Bas ¢ 
Cia € T Ger T Ss 3 dar ) 
cae ge une 


The generalized aerodynamic forces and moments can be ob 
The virtual 
work corresponding to a virtual displacement in any of the three 


tained through the use of the virtual work principle. 


types of motion is 


1 
6W; = b f ; II (h, a, &, 052) dx, 7 = 1,2,3 (11) 
where 
62, = 6h 
622 = b[x + (1 2) |ba (12) 
623 = (1 — x?)bt 


Substituting 10 and 12in 11 and performing the necessary integra 


tions yields: 











_ f h f 
bW, = —pV2bbh | 21 ao» + doa @ + Gog (r) + 
b b 
, h c 7 
4t ay Fi + Ug a + aig h 2 
h : 
iW, = p10 (au + dia @ + aig EG 
h “bh 4 
" (13) 
of. a + as - ( : ") 
la $ h ' } 
a ) h j 
sW; = p V2b8e | 2( a b + doa a + dog ;) (£) + 


( he ' ¢ 3n 4 ¢ 1 
a tT Qig ¢ t 3c — 
ss b ~—s ™ b 8 . 3 a b bad 


If now the right-hand sides of Eqs. 13 are multiplied and divided 
by b%w?, if the time-dependency factor is reintroduced, and if the 
definitions of equations (7), (8), and (9) are inserted, the following 
expressions are obtained for the generalized forces: 


6W, — h ¢ 
L= ah = rpbiwre' (1, b +Laa+Le ‘) 
bW, = h e 
M = “a = mpbtw%e (a, b + Maa + M; 4 (14) 
: 6W3 ees Pe 7 a 
N = a = rpb3we'! (a 5 + Naa + Ne 4 
where 
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Fic. 1 
L, = 1 i- C(k 
k 
l | a 2 
ba * ! ‘- Ce C(k 
2 k k k 
> l » 
Ly = i- C(k) 4 C(k 
{ k k 
VU, = 1/2 
Ma = (3/8) i(1/k (15) 
3 
We —_ ; ¢ t 
8 2k k 
} l 
N, = 1 C(k 
{ k 
; i. i 
Ne = i t= Ce C(k) 
8 k k kA 
7 = l 
Nr = i C(k) 4 + C(k) 
12 2k 2k? k? 


and where / is the oscillatory lift on the airfoil, 1/7 is the oscilla 
tory moment about the airfoil quarter-chord, and .V corresponds 
to a second moment about the airfoil mid-chord point. In the 
(14) 


The coefficients 


expressions for the forces and moments of Eqs ind (15), 
the notation of reference 2 has been adhered to 
La, Lea, My, and Mg are tabulated as functions of 1/k in reference 
2. The remaining coefficients Ly, Wr, Ny, Na, and N¢ are tabu- 
lated in reference 3 

In order to introduce the camber change coefficients in a prac 
tical flutter analysis, it would probably be most convenient to use 
the normal modes of vibration of the wing as degrees of freedom 
However, the chordwise, as well as the spanwise, distribution of 
amplitude in each normal mode would have to be known. In 
each mode and at each of several stations along the wing span, the 
actual deformation of the mean chord line would be plotted and 
then reduced to equivalent h, a, and ¢ components. The aero 
dynamic terms would thus be treated by strip theory, although 
the mechanical terms could be handled on a ‘‘three-dimensional”’ 
basis. Such a flutter analysis procedure might be highly ap 
proximate, but it would represent a relatively simple initial ap 
proach to the problem of accounting for the effects of camber 
changes in flutter. 

A rough indication of the magnitude of the lift and pitching 
moment due to parabolic camber change can be given by replac 
ing the parabola with a triangle, as shown in Fig. | This tri 
angular deflection pattern corresponds to a wing that translates 
an amount h, rotates through an angle a about the quarter-chord 
point, and has a 50 per cent chord control surface, hinged at its 
leading edge, which deflects through an angle 8 relative to the 
In order that the maximum deflections in the two cases 


wing. 
may be the same, it is necessary to set h/b = '/.,a = 1,8 = 2, 
and ¢/6 = 1. The oscillatory lift and moment per unit span in 


the two cases are, then: 


fg 
N 
| 


r pb3w%e' [L(h/b) + Laa + Lg] 


: (triangular 
rpbtw? &@[M(h/b) + Maa + Mp8] sil 


= 
ll 
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L = rpbw* ew Lt ¢/b)) zs 
M = rpb‘o?  [M(¢/b)] (parabolic) 


For purposes of comparison, one need examine only the terms in 
the brackets in each case—i.e. 


(1/2)L, + La — 2Lg \,,. , 
(1/2)M, + Ma — 2Mgf “rianeular 


and 


Ly | 
Mz) (parabolic ) 


Using the formulas of reference 2 for Lg and /g and Eq. (15) for 
the other aerodynamic coefficients, one obtains 


1 
Li + Le — 2Lg = 0.576 — 0.72745 C(k) + 


l 
1.273 zp CO) (16) 
I 0.75 - C(k) +2 . C(k 
«= 0.40 — 1 e) + 2 (Rk) | 
k k2 | 
l oe: 1 
M, + Me 2Mg = 0.288 + 0.4247 + 0.637 
2 k k? 17) 
7 y , (1¢) 
My = 0.875 + 0.57 - 4 
k k2 


The two expressions for the lift and moment are of identical form, 
differing only in the numerical coefticients of each term. 

In case a more accurate realization of a particular camber de- 
formation pattern is desired than can be obtained with either the 
parabolic or two-segment approximation, two alternatives are 
available: (1) Derive the aerodynamic terms for higher order 
modes of camber deflection using the general method presented 
above; or (2) Employ more straight-line segments. Both possi- 
bilities have the drawback of adding considerably to the compu- 
tational labor involved in performing a flutter analysis. 
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Note on the Use of Two-Dimensional 
Compressible Flow Aerodynamic Coefficients 
in the Flutter Analysis of a Tapered Wing 


R. J. Werdes 
Design Engineer, McDonnell Aircraft Corporation, St. Louis 
March 20, 1953 


F appoman IS OUTLINED which makes it possible to compile, in 
a compact table for use with little or no subsequent inter- 
polation, the various compressible flow aerodynamic coefficients 
used in the flutter analysis of a tapered wing 

In the practical analysis of the flutter characteristics of a 
tapered wing, the span is usually idealized as a composite of con 
stant chord sections. This or similar treatment is necessary when 
considering the effect of spanwise variation in elastic deformation 
and the effect of surface taper. The two-dimensional aerody 
namic coefficients defining the forces and moments at these 
sections are functions of the reduced wave length 1/k defined as 
v/bw. where v is the forward velocity, b is the section semichord, 
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TABLE 1 
Values of 1/k for R = 2'/s 


0. 50000 2.0000 8.0000 
0.54525 2.1810 8.7241 
0.59460 2.3784 9.5137 
0.64842 2.5937 10.375 
0.70711 2.8284 11.314 
0.77111 3.0844 12.338 
0.84090 3.3636 13.454 
0.91700 3.6680 14.672 
1.0000 4.0000 16.000 
1.0905 4.3620 17.448 
1.1892 +. 7568 19.027 
1.2968 5.1874 20.749 
1.4142 5. 6569 22.627 
1.5422 6.1688 24.675 
1.6818 6.7272 26.909 
1.8340 7.3360 29 344 


32.000 


and w is the circular frequency of oscillation. Thus, for a value 
of v/w constant over the entire wing, a different value of 1/& is 
defined at each section. While almost any convenient sectioning 
of the span is adequate when considering the effect of elastic de 
formation, a judicious sectioning of the span will greatly simplify 
the calculations when considering the effect of surface taper 
Suppose the two-dimensional aerodynamic coefficients have 
been tabulated corresponding to values of 1/k which form a 
geometric series of ratio R. If the tapered wing under consider 
ation is then idealized as a composite of constant chord sections 
oriented spanwise such that the section average chords also form 
a geometric series of ratio R, the aerodynamic coefficients re- 
quired in the evaluation of the oscillatory forces and moments 
can be obtained directly from the prepared table. In particular, 
for R = 2 
for the values of 1/k given in Table 1 
tively few values of 1/% are required in covering the range of usual 
Fig. 1 represents the tip portion of a nondimension 


’, the aerodynamic coefficients would be tabulated 
It will be noted that rela 


interest. 
alized tapered wing of infinite span, the chords of which are nor 
malized with respect to the tip chord, and on which a standard 
For R = 2, the 
1 with eight sections contained in 


span is defined for a taper ratio of one-half. 
sectioning is as shown in Fig 
each portion of span having a taper ratio of one-half. For a given 
taper ratio, the number of sections depends on the spanwise 
In most cases the root 


How 


ever, the aerodynamic coefficients applicable to this partial sec 


location of the normalized root chord 
chord will not coincide with a standard section boundary 


tion can be conveniently written in terms of the adjacent section 
coefficients. The following relations are used in locating sections 


on the given wing span: 
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infinite span showing spanwise sectioning scheme. 
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"1 = {94 C1) — 1 
3 = R-1-1(4 = 1,2,3. 


where 
» = span reference dimension 
R = ratio of adjacent section chords 
C = chord 
= subscript denoting section boundary location with 7 = 
l at the wing tip 
r,t = subscripts denoting wing root and wing tip locations, 
respectively 
For R = 2%, the values of »; can be obtained from Table 1 by 
subtracting unity from the values of 1/& starting with 1/k = 1 
When the value of (», — »,) has been established, all dimensions 


can be converted to full-scale wing dimensions 
The scheme outlined above is ideally suited to punched-card 
computational methods, since the functions can be stored with 


relatively few cards and practically no interpolation is required. 


An Observation on the Vortex System of Dual- 
Rotation Propellers 


J. Eisenhuth and B. \W. McCormick 
Research Associates, Ordnance Research Laboratory 
The Pennsylvania State College, State College, Pa 


March 5, 1953 

I‘ ORDER TO DETERMINE the induced effects associated with a 
pair of dual-rotation propellers, it is commonly assumed 

that the vortex system aft of the propellers is formed by the 


intersection of right- and left-handed nondistorted helical vor 


tex sheets.' Recent water-tunnel* experiments have shown 
that the above concept might be in considerable error. The 
accompanying high-speed photograph, Fig. 1, shows a dual 
rotation propeller operating on the tail of a torpedo model. The 


outer edge of the vortex sheet is easily discernible because of the 


occurrence of cavitation in this low-pressure region. In addition 


to the use of high-speed photographs, the behavior of the tip vor- 


* These tests were conducted in the 48-in. water tunnel of the Ordnance 


Research Laboratory, operated under Contact NOrd 7958 with the Bureau 


f Ordnance, Department of the Navy, by The Pennsylvania State College 
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tices was viewed with stroboscopic light. These photographic 
and visual observations are extremely enlightening with regard 
to the qualitative nature of the combined trailing vortex system. 
The tip vortices trailing from the individual blades of each pro- 
peiler intersect in the plane of the after propeller. The flow 
configuration is apparently unstable at these points of intersec 
tion. Instead of the vortices from each propeller blade main 
taining their expected helical shapes, the ones from the forward 
propeller break and join with those from the after-propeller, 
series of closed loops that move downstream with 


forming a 
In view of the observed behavior of the tip 


periodic spacing. 
vortices, it is reasonable to assume that the inner portion of the 
vortex sheets must also undergo some similar distortion, thus 
invalidating the assumption of true helical surfaces 
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A Note on Second-Order Supersonic Flow 
Theory 


Robert E. Lavender 
Technical Feasibility Studies Office, Ordnance Missile Laboratories, 
Redstone Arsenal, Huntsville, Ala. 
January 21, 1953 
fo DYKE HAS RECENTLY POINTED OUT that the solution to 
the second-order equation for the supersonic flow past non- 
lifting bodies of revolution whose slope of the meridian line has 
discontinuities requires a superposition of five basic solutions to 
the first-order equation.!' The writer feels that a discussion of 
the suitable strengths of these solutions is warranted. 
Consider a body of revolution (Fig. 1) which has discontinuities 
in curvature or slope at points VY, to Y The first-order solu 


tion can be written as follows: 


n—l 4 
gn = —AX,b(t,) — >> BAX, — &)?al Ti") - » ( 
i=0 y l 
4 
(XX, — Eny)/? g(Ts ~ > Diy VXn- h(t (1) 
y l 
where C;., = Diy = 0 for Xn < Xiy 
, l 
¢r, = ABelt + 8 BAX, — &) c(Ti") 4 
1=0 
1 4 
8B >. ¢ VX, — Exy UT )+~p >> D 
oak y=1 
(Xn — Eny) > Rk ) (2 
n l 
or, = —Ad(tn) — DS BAX. — &)(Ti") —- 
1=0 
- £.. 
> ( vz -_ , hi 5 )— ke D 
1 y=1 
(XY, — Emy)~ “277 (T (3) 
-1 
vr = —A(Y,) (t m 3 Bd(T;") — 
’ 7=0 
4 
YD, Ca Xn — Sey)" UT 
y=1 
& 
> DilXi — bu KT (4) 
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Fic. 1. Body of revolution. 


As noted by Van Dyke,! the other second derivatives can be ex- 
pressed by: 


Grrn = Bgrrn — (Grn/Rn) (9) 


Grr, = R,” — Ra’ Grr, (for o¢,, = R;’) (6) 


Eqs. (1) through (6) give the first-order potential and its deriv- 
atives at station Y, along the body, where 


t, = BR,/X, 


te = BR’, Xo = 1, Ro = Ro’, Ro” = 0 
T;" = BR»A/(Xn — &:) 
Si = Xi = BR; 


The functions a(t) through m(¢t) are defined and tabulated in the 


referenced report. 


Taking the first-order tangency condition to be ¢g,, = R,,’ (as 
suggested by Van Dyke!), 


A = ¢?/V 1 — Bre”, « R,’ (7) 


1 
- | R,’ — ABe(t,) — 
BXn — En-1)e(T n-1") 


n—2 4 
/ 
B > BAX, — ETI") — BOD Coy VXun — Exy (Tey) — 
y=1 


1=0 


5B, ; = 


1 
B > DilX, = ty) a) (8) 
y=1 


Diy = (Riyo’ — Riya’) * BR B (9) 
Vv BRxy 
Cry = Riyo” — (¢er)kua — Reyo'(¢rr)kva + 
sB- B*Riyp’ ? ' y 
(Riyr’ — Riya’) [7 3 
“ Riy’ — (10) 
Ris 8 88 
where Ry,q" = R” just ahead of the discontinuity at station V,, 


and Rx,” = R” just behind the discontinuity at station Y;,,. If 
the discontinuity at station Y,;, is in curvature only and not in 
slope, Eqs. (9) and (10) reduce to 


Diy = 0 
Vv BR ky 


Cie = = . Chins” _ Ria 3 
B = B?R ky 


Eqs. (7) through (10) are the constants defining the strengths of 
the basic solutions which are needed to obtain the first-order 
potential and its derivatives. 


A particular solution, y, to the second-order perturbation 
equation has been obtained by Van Dyke in terms of the first- 
order potential, yg, and its derivatives. The second-order po- 
tential, ¢, consists of this particular solution plus a complemen- 
tary solution, x, to the first-order equation, so that (as given by 
Van Dyke’) 


Xrn = Ril + orn) — Ven (11) 
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The expressions for x; and x, can be written: 


n—1 


Xrn = EBe(tn) + B DY Fi(Xn — &)c(TH") + 
1=0 
4 4 
BY) GilXn — Envy) '7R(Tiy") +B > Hil X,- 
y=1 y=1 


4 ‘ fepieacieccms 
Eey)~"*m(Tay") + DS B Sky VXn — Ey Tay") (12 


y=1 
where Gy, = Hiy = Jiy = O for Xn < Xuy and 
n—1 
Xen = —Ed(tn) — YS F(Xn — &)b(Ti*) — 
1=0 


4 4 
> » Gil Xn — Eey)~'/9( Tey") + Zz. Hu Xn — Eay)~'/ UT yy")- 
y=1 y=1 


4 oi . 
} Jey VXn — Sey (Tey) (13 
y=l 
Eqs. (12) and (13), along with the values of y,, and y,, (ob- 
tainable! in terms of the first-order values), define the second- 
order perturbation velocities at station X, along the body: 


Prn - Vrn ss Xrn (14 
orn = brn + Xzn) 
Taking the second-order tangency condition to be given by 
Eq. (11), 


E= [Ro’(1 + ¢2,) — ¥ro] Be( ty) (15 


| eva + ¢r,) 


F,-1 = - 
MXa = ba-1)0(T a=1*) 
n—2 
Vin — EBe(tn) — B D>. F(X, — E)e(TI") — 
1=0 
4 4 
BD) Gal Xn — Fey)“ '4R(Tay") — B DD Hel Xn — gxy)~”? X 
y=1 you 


4 
m(Tky") — B >> Siy VXn — Ety ite) (16 


yv=1 


In determining the constants defined by Eqs. (8) and (16) when 
X, corresponds to a station Y,, with a discontinuity in R” or R’, 
take X, = Xiya. 

hp = (Wyn — Wkva) VBRiy (17 


! V BR; ; , Reyr’ 
Gy = . (Riyo’ — Reva (1 T Priya — sid 
B B 


3 
(Wiyo — Wkva) (18 
8Riy , 


Fix = Gu (Cu Dy), D,. QV 
= () . Beg = B 


(Urey vd Vrkya) wi 


(19 


If the discontinuity at station Y,, is in curvature only and not 
in slope, Eqs. (17) and (18) reduce to 


Ay = 0 
Gry = (—-VBRiy/B) (Yr vp — Vrava 


Eqs. (15) through (19) are the constants defining the strengths of 
the basic solutions which are needed to obtain the complemen- 
tary solutions from Eqs. (12) and (13). The generalization of the 
above equations for bodies with more than four discontinuities 
is easily seen. 

As pointed out by Van Dyke, the proper expression for the 
second-order pressure coefficient is given by the exact expression: 


2 ~1 
Con = (4 re M? x 
yir7\) 2 


[1 — (1+ eq)? — den?) 6 b-1 


where the second-order perturbation velocities are given by Eqs. 


(14). 
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Van Dyke 


influence Matrix of Beam Vibration 


L.$.Dzung 
Brown, Boveri & Cie., Baden, Switzerland 
February 13, 1953 

HE DIFFERENTIAL EQUATION of beam vibration (EJ y")” = 
eae (E = Young’s modulus, / = area, moment of inertia, 
y = deflection, « = mass per unit length, w = angular frequency, 
and primes denote differentiation with respect to beam length) 
when applied to a uniform beam section (J and yu constant) is 
readily solved by standard methods and leads to the following 
matrix equation :1 


Vi+1 j a b € vi 

0541 _ uc 7 a b : 0; 
Mi+1 ub uc fj a M; 

Seva | ua ub ue FZ | Ss 


where 0 = y’, M = Ely’, S Ely’’’ denote the slope, the bend 

ing moment, and the shearing force, respectively. 

andi + 1 denote the left and right end of the section 
The elements of the influence matrix are: 


Subscripts 7 


7 = (1/2) (cosh a + cos a) 
a = Ela = (1/2) (sinh a + sin a) (1/a) 
b = Elb = (1/72) (cosh a — cos a) (1/a)? 
c = Elé = (1/2) (sinh a — sin a) (1/a)3 
“= pw* 
a= pw?/ El 
If the section length / is chosen short enough so that at << 1, 
the elements can be approximated by 
j = ]°/O! = 1 
a= Efé =['/1! =] io 
b = Elb = 22/2! = 2/2 ad 
c = Efe = [*/3! 13/6 


This is the result obtained by Morton.? 

If the section is not uniform but if the changes of 7 and yu are 
small and can be approximated by straight-line variations along 
the beam length, then the differential equation can still be solved 
by series integration. The result is similar to Eq. (1), except 
that it is more convenient to use the quantities y’ = W = 


W/E, y'"" = 5 = S/EI 


Eq. (1) becomes 


Via J a b ( Vi 
9:4 ! _| j a b 6: (3) 
Mis: ub ue j a M; 

| Si41 ud ub ue fF §, 


This equation is more symmetrical than Eq. (1) and, within the 


» 


approximation of Eqs. (2), is applicable to a nonuniform beam 





he 


mn: 


as well, provided 


AI/I = 2 Ap/p<<k 


where A denotes the difference between the values at station 7 
and at station 7 + 1. 


2(4 + n)/(2 — n) 


k = 


l or 2—i.e.,k = 10 0r ~. It is understood that the val- 
ues of J and yu just to the right of station 7 are taken to calculate 
the matrix elements, while the values at the respective stations 


are taken for Wand S. 


form = 
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If a concentrated mass m is attached to the left end just to the 
right of station i, then the shearing force of the whole section is 
increased by mw*y. The first column of the matrix in Eq. (1) or 
(3) must be increased by mw? multiplied by the respective last 
column. Fora massless beam, u = 0, this is the system of equa 


tions used by Prohl,* Myklestad,* Targoff,> and Thomson.® 


If the same mass is attached to the right end just to the left 
of station 7 + 1, then the shearing force at station i + 1 is in 
creased by mw*y;+;. The last row of the matrix must be in 
creased by mw? multiplied by its first row. For a massless beam, 


this is the equation given by Thomson. ! 


A word of caution is in order if, for the calculations of natural 
frequencies and natural modes, the entire beam is divided into 
sections to each of which Eq. (1) or (3) is applied. If the cross 
section of the beam varies smoothly—e.g., in case of an airplane 
wing—then 1/ and S or WM and S are usually continuous at the 
joint between two sections, except due to concentrated mass or 
external constraint. If the beam cross section changes abruptly at 
e.g. in case of a rotating shaft with stepwise vary 


then J/ and S are still continuous usually but not 


a Station point 
ing diameters 


M and S due to abrupt change of / 
propriate to use Eq. (1) than Eq. (3), in order that the right end 


In this case it is more ap 


conditions of one section can be used directly as the left end con 


ditions of next section. 
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Regarding the Computation of Unsteady Air 
Forces by Means of Mathieu Functions 


Henry E. Fettis 
Mathematician, Flight Research Laboratory, Wright Air 
Development Center, Wright-Patterson Air Force Base, Ohio 


February 13, 1953 


o- THE PUBLICATION OF the last letter on this subject by 
the author,! significant new results have been obtained from 
calculations performed by Dr. Gertrude Blanch, at the Institute 
for Numerical Analysis of the National Bureau of Standards at 
Los Angeles, based on the equations of Reissner*® and Haskind.‘ 
Some preliminary results of these computations are shown in 
Table 1, together with the corresponding values by Dietze’s 
method.? From these results it can be seen that essentially the 
same values are obtained by the exact solution in terms of 
Mathieu functions as were obtained by the approximate method, 
the effect of the approximation being confined to discrepancies 
(in general) in only the third significant figure. Thus it may 
be safely concluded both that Dietze’s results are fundamentally 
correct and that the method of Mathieu functions can be made 
to duplicate results previously obtained by other means to a satis 
factory degree of accuracy. It is also noted that the results ob 
tained from the equations of Haskind and Reissner satisfy the 


relation 
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TABLE | 





Fp iota ~ Ly ——_——- 
w a b 
0.6 —0.05032 —2. 289197 — 0.0550 —2 . 29227 
LT +0.03390 — 1.993852 +0.02715—1.9977% 
aa —- Le — 
w a b 
0.6 — 4.64460 — 1.80363: —4.6581—1.79677 


0.7 


—3.55325—1. 595231 


La + M, = [1 — (t/w) JL, 
identically with the exception of one unit in the fifth decimal place 


for w = 0.6, as shown in Table 2. 


TABLE 2 


w La + Mh {1 = (1 w)| La 
0.6 —3.86564—2.20534: —3.86564 —2. 20533: 
( — 2.81446 —2.04228: 


).7 —2.81446—2.04228:1 


In conclusion, the author wishes to say that the pointing out of 
the discrepancies that exist between the values of the coefficients 
as presented by the Messrs. Timman and van de Vooren and 
those otherwise computed is in no way intended to discredit 
either the validity of their method or the value of their work. It 


On the Measurement of Surface Pressure with 
a Static Probe 


W. W. Willmarth 
Graduate Student, Guggenheim Aeronautical Laboratory, 
California Institute of Technology, Pasadena, Calif. 


March 2, 1953 


Sorrace PRESSURES on two-dimensional bodies of small thick 

ness ratio have been measured by traversing a conventional 
static pressure probe over the surface. The method can be ap- 
plied to thin bodies that are difficult, if not impossible, to instru- 
ment with conventional surface pressure taps. 
should prove useful in problems involving the shaping of aero- 
dynamic surfaces to obtain a given pressure distribution, for ex- 


The same method 


ample, in wind-tunnel tests on pylons, fuselage-wing fillets, etc. 


The fact that surface pressures on thin bodies can be measured 
with precision at points off the surface depends on two effects. 
The first is the effect of viscosity in providing a boundary layer 
through which the static pressure is constant in a direction normal 
to the body surface when the surface curvature is small. The 
second effect depends on the fluid being irrotational outside of the 
boundary layer. In conventional notation, using two-dimen 
sional perturbation velocities, « and v, the irrotationality con 
dition is 

Uy — tr = O (1) 
The linearized boundary condition on the body surface is 
v/U.,, = (body slope) (2) 
on y = 0, and the pressure coefficient is 
Cp = —(2u/U,,) (3) 


From Eqs. (1) and (2) it is apparent that in the inviscid flow 
near the surface the change in wu in the y-direction is proportional 


Comparison of Lift and Moment Coefficients for 1 = 0.7 (a) Values by Reissner’s Formulas; (b) Values by Dietze’s Method 


— 3.5630 —1. 58307 


Fic. 1. 


——_—_—__— — M, — —— 

a b 
0.778961 —0.401717 0.7806 —0. 4053: 
0.738792 —0.447051 0.7394 —0.4535i 

Ma. - ; 

a b 
— 0.42508 —2 . 851831 —0.4258—2.86171 
—0. 45912 —2. 43883: —0. 4691 —2. 4434: 


is, rather, to call attention to the presence of some (perhaps ob 
scure) numerical error that has somehow entered into their 
calculation scheme. 
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“SLIDING PROBES 


Sketch of the probe configuration for a 3 per cent thick 
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to the x-derivative of the slope of the body when the boundary- 
layer growth is small. For flat surfaces 0u/Oy is negligible, and 
from Eq. (3) the pressure change normal to the surface (y = 0) 
will be extremely small. On slightly curved surfaces, the change 
of body slope is small and the normal pressure gradients must also 
be small. 

Using the above considerations, a pair of slender static pressure 
probes aligned with the chord of a thin two-dimensional airfoil 
have been used to obtain surface pressure distributions at tran 
sonic speeds and small angles of attack. The probe arrangement 
on the airfoil cross section is shown in Fig. 1. In order to mini 
mize tunnel blockage, the probes slide inside larger supporting 
tubes, which lead back from the airfoil and out through the tunnel 
wall. The probes were held together by a loop of wire approxi 
mately two diameters from the tips. 
probe boundary layer turbulent, which is important when meas 
The ratio of probe diameter to 


The wire loop also made the 


uring steep pressure gradients.! 
airfoil thickness was of the order of one-half. 


The effect of the small three-dimensional probe disturbances 
on the basic two-dimensional flow over the airfoil was checked 
with a Mach-Zehnder interferometer.? An interferogram of the 
flow over one surface of a 8 per cent thick airfoil was made and 
the pressure distribution was computed. The probes were then 
installed and used to find the surface pressure distribution. The 
airfoil was investigated for an angle of attack of approximately 
| deg. at two supersonic Mach Numbers. In one case, the bow 
shock wave was attached with completely supersonic flow over 
the airfoil. The second case was at a lower Mach Number with 
a detached bow wave giving subsonic flow over the front portion 
of the airfoil. The agreement between interferometer and probe 
measurements was good except in the vicinity of the airfoil nose, 


as shown in Fig. 2. 
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Using this probe technique, pressure distributions on a 5 per 
cent thick double wedge have been measured at high subsonic 


speeds.* 


REFERENCES 
! Liepmann, H. W., Roshko, A., and Dhawan, S., On the Reflection of Shock 
Waves from Boundary Layers, NACA TN 2334, p. 20, April, 1951 
2? Ashkenas, H., and Bryson, A. E., Design and Performance of a Simple 
Interferometer for Wind Tunnel Measurements, Journal of the Aeronautical 
Sciences, Vol. 18, No. 2, p. 82, February, 1951 
> Cole, J. D., Solomon, G. E., and Willmarth, W. W., Transonic Flow Past 
Simple Bodies, presented at the Twenty-First Annual Meeting, IAS, New 


York, January 26-29, 1953 


Concerning ‘‘The Subsonic Calculation of 
Circulatory Spanwise Loadings for Oscillating 
Airfoils by Lifting-Line Techniques’’' 


M. A. Dengler and Martin Goland 
Midwest Research Institute, Kansas City, Mo. 
February 19, 1953 


—— OF OUR PAPER shows that we do not make entirely 
clear the nature of the mode shapes used for the swept-wing 
Our intent was to study the effects of geometrically 
Hence, after a calculation was performed 


calculations. 
sweeping a plan form. 
for an unswept wing, a similar computation was made for a 45° 
swept wing having the same geometric downwash along the rear 
quarter-chord line as the unswept case. The paper is clear con- 
cerning the geometric downwash for the unswept wings; the 
geometric downwash boundary conditions for the swept wings 
are taken as identical. 

We should also like to call attention to an inadvertent omission 
from our bibliography of the valuable series of papers by W. P. 
Among these are the Aero. Res. Council Tech. Repts., R. 
The last paper bears 


Jones. 
and M. Nos. 1920, 2142, 2215, and 2470. 
particularly on our approach and presents a method by which 
Falkner’s discrete-vortex scheme for calculating surface loadings 
can be made part of a technique for air-load calculations on os- 


cillating, finite-span airfoils. 
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Errata—‘‘The von Karman Integral Method as 
Applied to a Turbulent Boundary Layer’’’ 


Nguyen Van Le 
Preliminary Design Engineer, AiResearch Manufacturing 
Company, Los Angeles, Calif. 


Februery 23, 1953 


i THE EQUATION PRECEDING Eq. (5), in Eqs. (5) and (6), and 
in the last equation, there should be 


6 O 
—u'y’ dy 
0 Ox? 


instead of 
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Errata—'‘‘Elastic Stability of Simply Supported 
Flat Rectangular Plates...’’' 


Robert G. Noel 
— Laboratory, North American Aviation, Inc., Downey, 
Calif. 


March 18, 1953 


¥ igs AUTHOR IS GRATEFUL TO Dr. Gorden C. k. Yeh, Project 
Engineer, Analytical Engineering Department, Reed Re 
search, Inc., for pointing out that certain corrections should be 
made in reference 1. 

These typographical errors were made in the manuscript sub- 
mitted for publication and therefore are not the fault of the pub- 
lisher. It should be noted that the validity of the result curves 
given in the paper is not affected. 
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(1) In Eq. (3), the first term of the integral should read: 


J ( o*w n =) 
i\ ax? oy? 


(2) In Eq. (5), (0?w/dx?)? and (07w/Oy?)? should be replaced by 
(Ow/Ox)? and (Ow/Oy)?’, respectively. 

(3) In the first term of Eq. (6), the number 16 should be re. 
placed by the number 8. 

(4) In Eqs. (7) and (8), the quantity [1 — n8?]? should be 
replaced by the quantity [1 + n?8?]?. 

(5) The factor outside the bracket in Eq. (9) should read 
SkprB2at/2?A 

(6) Expression (10) should be set equal to zero. 

(7) In the expression for A following Eq. (10), the first term 
should read (1 + 778?)? 

(8) In the fourth line under the section ‘“‘Use of the Curves,” 
‘*Figs. 3’’ should be replaced by “‘Figs. 2.”’ 

(9) The expression for kp in Fig. 4 should be divided by §?, 

(10) In the sentence next to the last sentence in the paper, 
“*R,” should be replaced by “‘R;.”’ 
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